Robust Control of a Multivariable Experimental Four-Tank System

The goal is to explore some results of the two papers
regarding the analysis of the properties of a multivariable Four-Tank System with 2 inputs and 2 outputs. We will also provide controllers for robust stabilizing the linear tangent model of the non linear system around the origin. Our dissertation will follow these keypoints:
The goal is to control the levels and in the two lower tanks using the two pumps. The plants control inputs are the two input voltages of the pumps, and . The fòow produced by pump i corresponding to a voltage is . The i-th valve distributes the output ow of the pump by sending to the lower tank and the remaining to the upper tank. Clearly . The gravity acceleration is g while is the cross section area of tank i and is the cross-section area of the outlet hole of each tank. A sensor provides an output voltage proportional (with constant ) to the water height , therefore the two measured outputs (two sensors with same ) are
Denoting with i the height at the equilibrium in each tank, the time costants are given by:
The values of these parameters are taken from the Johansson' paper.

1. Poles and Zeros location and Minimum Phase property analysis.

The control object of our analysis is the MIMO transfer futnction matrix of the linearized model around an equilibrium point:
= with and
Let us first of all compute the poles and the zeros.
We know that in a MIMO system the tranfer function's poles can be computed as the roots of the pole polynomial. The pole polynomial , i.e. the least common denominator of all non-zero minors of all orders of G(s).
The minors of order 1 are the entries of the matrix, while the only minor of order 2 is the transfer function determinant:
.
so .
Then we can state that the poles of the transfer function are each of them with multeplicity equal to 1.
Similar considerations can be done for the computation of the zero. In the general case, the zeros can be computed as the roots of the zero polynomial i.e. the greatest common divisor of all the numerators of all order r (normal rank) minors of G(s), provided that these have been adjusted so to have the pole polynomial as denominator. So if the normal rank of is 2, we need to check only the minor of order 2, i. e. the determinant of the matrix.
Additionally, becuse of the fact that the transfer function is square (p = q = 2), we could directly say that the solutions of 0 gives the transmission zeros which do not coincide with a pole.
Let's check when the determinant is different from zero. If it is equal to zero only at few values (normal rank = 2), we can state that such values are the transmission zeros of the transfer function.
so:
a polynomial of order 2 whoose roots are:
, with , a = ,
and consequently the term = ( - )
So, for the physics of the system, parameters will never be such that these two solutions are complex, and they costitutes the only two values with respect to which the matrix G(s) loses rank. So G(s) has normal rank equal to 2, and that solutions are the symbolical two trasmission zeros of the trasfer function.
Let us validate our computations using the calculator:
clear all
close all
clf
 
syms g1 g2 k1 k2 n1 n2 t1 t2 t3 t4 s A1 A2 A3 A4 kc % inizializing the symbolic parameters of the system
 
G = [(g1*k1*n1/(1+t1*s)), ((1-g2)*k2*n1)/ ((1+t3*s)*(1+t1*s)); ((1-g1)*k1*n2)/ ((1+t2*s)*(1+t4*s)), (g2*k2*n2)/(1+t2*s)] % inizializing the plant
G = 
d = det(G)
d = 
eqn = d == 0;
zi = solve(eqn, s, 'ReturnConditions', true);
simplify(zi.s)
ans = 
These are the result we obtained, barring semplifications.
From the above analysis, it is easy to verify the properties of those zeros.
Thinking about the physical interpetation these results make sense. We have a non minimum phase system if , i.e. if the entering flow into the two lower thank ( the ones with respect to which we want to regulate the water level) is smaller than the total one entering the two upper tanks. Viceversa, we have minimum phase property is the upper tanks receive more flow tha the others.
If instead the total flow going to the left tanks 1 and 3 is equal to the total flow going to the right tanks 2 and 4, we have to manage with some more difficulties derived by the presence of the zero in 0.
Note that:
.
So, in the case in which and , we have that and , so the two zeros of the transfer function are close to be coincident with the two poles , .
Note also that if one of the two or both we have not transimition zeros in the transfer function.

2. Directions of the zeros and the poles, blocking property of the zeros.

In order to compute the zeros directions, we need a minimal realization {A,B,C,D}: without it, we could not use the useful results about the system matrix.
Let us check if the proposed realization in the Johansson's paper is minimal.
To do this, we need to check the observability and the controllability properties of the natural modes of that system. If the system is fully controllable and fully observable, we can state that the realizazion {A,B,C,D} is minimal.
So let us compute the Observability and the Controllability Matrices and check if they are full rank.
A = [ -1/t1, 0 , A3/(A1*t3), 0; 0 , -1/t2, 0 A4/(A2*t4); 0, 0, -1/t3, 0; 0,0,0, -1/t4]; % matrix A
C = [kc, 0,0,0; 0, kc,0,0]; % matrix C
B = [(g1*k1)/A1, 0; 0, (g2*k2)/A2; 0, (1-g2)*k2/A3; (1-g1)*k1/A4, 0]; % matrix B
D = zeros(2,2);
 
A_2 = A*A;
A_3 = A*A*A;
Obv = [C; C*A; C*A_2; C*A_3]; % observability matrix
rank(Obv)
ans = 4
Con = [B , A*B, A_2*B, A_3*B]; % controllability matrix
rank(Con)
ans = 4
Both matrices have rank = 4, equal to the dimention of the state. This implies that the system is fully controllable and observable. The state space representation is then a minimal realization of the transfer function G(s).
Let us now substitute the symbolical quantities with the values used in the Johansson paper and try to compute these quantities.
Our considerations on the zeros and poles directions will regard the operating points and and around which the system exhibits respectively minimum-phase (MP) and non-minimum-phase (NMP) characteristics.
A1_n = 28; % cm^2
A3_n = 28; % cm^2
A2_n = 32; % cm^2
A4_n = 32; % cm^2
c1_n = 0.071; % cm^2
c2_n = 0.057; % cm^2
c3_n = 0.071; % cm^2
c4_n = 0.057; % cm^2
kc_n= 0.5; % V/cm
g_n = 981; %cm/s^2
 
 
% MP parameters
 
h01_mp = 12.4; % cm
h02_mp = 12.7; % cm
h03_mp = 1.8; % cm
h04_mp = 1.4; % cm
 
k1_mp = 3.33; % cm^3/Vs
k2_mp = 3.35; % cm^3/Vs
 
g1_mp = 0.7;
g2_mp = 0.6;
 
t1_mp = (A1_n/c1_n) * sqrt(2*h01_mp/g_n); % s
t2_mp = (A2_n/c2_n) * sqrt(2*h02_mp/g_n); % s
t3_mp = (A3_n/c3_n) * sqrt(2*h03_mp/g_n); % s
t4_mp = (A4_n/c4_n) * sqrt(2*h04_mp/g_n); % s
 
n1_mp = (kc_n* t1_mp)/A1_n;
n2_mp = (kc_n* t2_mp)/A2_n;
 
% NMP parameters
 
h01_nmp = 12.6; % cm
h02_nmp = 13; % cm
h03_nmp = 4.8; % cm
h04_nmp = 4.9; % cm
 
k1_nmp = 3.14; % cm^3/Vs
k2_nmp = 3.29; % cm^3/Vs
 
g1_nmp = 0.43;
g2_nmp = 0.34;
 
t1_nmp = (A1_n/c1_n) * sqrt(2*h01_nmp/g_n); % s
t2_nmp = (A2_n/c2_n) * sqrt(2*h02_nmp/g_n); % s
t3_nmp = (A3_n/c3_n) * sqrt(2*h03_nmp/g_n); % s
t4_nmp = (A4_n/c4_n) * sqrt(2*h04_nmp/g_n); % s
 
n1_nmp = (kc_n* t1_nmp)/A1_n;
n2_nmp = (kc_n* t2_nmp)/A2_n;
 
 
% MP state space inizialization
 
A_mp = eval(subs(A, [A1, A2, A3,A4, t1, t2, t3 ,t4], [A1_n, A2_n, A3_n, A4_n, t1_mp, t2_mp, t3_mp, t4_mp]));
B_mp = eval(subs(B, [A1, A2, A3,A4, g1, g2, k1 , k2], [A1_n, A2_n, A3_n, A4_n, g1_mp, g2_mp, k1_mp, k2_mp]));
C_mp = eval(subs(C, kc, kc_n));
D_mp = eye(2)*0;
 
G_sys_mp = ss(A_mp, B_mp, C_mp, D_mp);
 
% NMP state space inizialization
 
A_nmp = eval(subs(A, [A1, A2, A3,A4, t1, t2, t3 ,t4], [A1_n, A2_n, A3_n, A4_n, t1_nmp, t2_nmp, t3_nmp, t4_nmp]));
B_nmp = eval(subs(B, [A1, A2, A3,A4, g1, g2, k1 , k2], [A1_n, A2_n, A3_n, A4_n, g1_nmp, g2_nmp, k1_nmp, k2_nmp]));
C_nmp = eval(subs(C, kc, kc_n));
D_nmp = eye(2)*0;
 
G_sys_nmp = ss(A_nmp, B_nmp, C_nmp, D_nmp);
 
s = tf('s');
 
% inizializing the transfer function of the plant with MP property
 
Gs_mp = [(g1_mp*k1_mp*n1_mp/(1+t1_mp*s)), ((1-g2_mp)*k2_mp*n1_mp)/ ((1+t3_mp*s)*(1+t1_mp*s)); ((1-g1_mp)*k1_mp*n2_mp)/ ((1+t2_mp*s)*(1+t4_mp*s)), (g2_mp*k2_mp*n2_mp)/(1+t2_mp*s)]
Gs_mp = From input 1 to output... 2.61 1: ---------- 62.7 s + 1 1.41 2: ---------------------- 2709 s^2 + 120.3 s + 1 From input 2 to output... 1.5 1: ---------------------- 1498 s^2 + 86.59 s + 1 2.837 2: ----------- 90.34 s + 1 Continuous-time transfer function. Model Properties
% inizializing the transfer function of the plant with NMP property
 
Gs_nmp = [(g1_nmp*k1_nmp*n1_nmp/(1+t1_nmp*s)), ((1-g2_nmp)*k2_nmp*n1_nmp)/ ((1+t3_nmp*s)*(1+t1_nmp*s)); ((1-g1_nmp)*k1_nmp*n2_nmp)/ ((1+t2_nmp*s)*(1+t4_nmp*s)), (g2_nmp*k2_nmp*n2_nmp)/(1+t2_nmp*s)]
Gs_nmp = From input 1 to output... 1.524 1: ----------- 63.21 s + 1 2.556 2: ---------------------- 5128 s^2 + 147.5 s + 1 From input 2 to output... 2.451 1: ---------------------- 2466 s^2 + 102.2 s + 1 1.597 2: ---------- 91.4 s + 1 Continuous-time transfer function. Model Properties
Let's first of all compute the zeros of the system in the way explained in point 1, both for MP and NMP cases, i.d. seeing the values for wihich
det_mp = (k1_mp*k2_mp*n1_mp*n2_mp*(g1_mp + g2_mp + g1_mp*g2_mp*s*t3_mp + g1_mp*g2_mp*s*t4_mp + g1_mp*g2_mp*s^2*t3_mp*t4_mp - 1))/((s*t1_mp + 1)*(s*t2_mp + 1)*(s*t3_mp + 1)*(s*t4_mp + 1));
zi_mp = zero(det_mp) % computation of the 2 minimum phase transmission zeros
zi_mp = 2×1
-0.0580 -0.0172
det_nmp = (k1_nmp*k2_nmp*n1_nmp*n2_nmp*(g1_nmp + g2_nmp + g1_nmp*g2_nmp*s*t3_nmp + g1_nmp*g2_nmp*s*t4_nmp + g1_nmp*g2_nmp*s^2*t3_nmp*t4_nmp - 1))/((s*t1_nmp + 1)*(s*t2_nmp + 1)*(s*t3_nmp + 1)*(s*t4_nmp + 1));
zi_nmp = zero(det_nmp) % computation of the 1 minimum phase and the 1 non minimum phase transmission zeros
zi_nmp = 2×1
-0.0562 0.0128
As we can see, the signs of the zeros are coherent with what we expected. Let's now compute the system matrix and see if its zeros coincides with the zero of the transfer function.
We have to remember that since the system is square (p = q), left and right zeros coincides. Furthermore, since the realization is minimal, the zeros of the system matrix with respect to which it drops rank are ONLY the transmition zeros, there is no presence of invariant zeros.
syms s
sys_matrix = [s*eye(4)- A, -B; C, D] % inzializing the system matrix
sys_matrix = 
sys_matrix_mp = eval([s*eye(4)-A_mp, -B_mp; C_mp, D_mp]); % inzializing the system matrix MP
sys_matrix_nmp = eval([s*eye(4)-A_nmp, -B_nmp; C_nmp, D_nmp]); % inzializing the system matrix NMP
 
det_sys_mat_mp = eval(det(sys_matrix_mp)); % determinant of MP system matrix
det_sys_mat_nmp = eval(det(sys_matrix_nmp)); % determinant of NMP system matrix
 
eq5 = det_sys_mat_mp == 0; % finding values with respect to which the system matrix drops rank
eq6 = det_sys_mat_nmp == 0;
 
zi_mp= eval(solve(eq5,s)) % 2 MP zeros
zi_mp = 2×1
-0.0580 -0.0172
zi_nmp= eval(solve(eq6,s)) % 1 MP 1 NMP zeros
zi_nmp = 2×1
-0.0562 0.0128
Which proves our expectations. Let us further validate our results using the MATLAB function tzero() to see if there is matching of computations.
tzero(G_sys_mp)
ans = 2×1
-0.0580 -0.0172
tzero(G_sys_nmp)
ans = 2×1
-0.0562 0.0128
Verified that the computed two pair of zeros are right, we can now find the zero directions solving the two problems for every zero:
, with and
where is the system matrix computed at one possible zero, , , respectively the right and left generalized eigenvector associated to z.
The vector costitutes the zero state direction , the zero input direction, the zero output direction.
We can reformulate this two problems as a generalized eigenvalue problem:
with and , and find easily.
Let us follow this approach.
M = [eye(4), B_mp*0;C_mp*0, D_mp]; % inizialization of M
 
G_mp = [A_mp,B_mp; C_mp, D_mp]; % inizialization of G (MP zeros)
G_nmp = [A_nmp,B_nmp;C_nmp, D_nmp]; % inizialization of G (1 NMP zeros)
 
[na,za,ma] = eig(G_mp,M) % computing n and m
na = 6×6
0 0 0 0 0.0000 0 0 0 0 0 0 0.0000 0 1.0000 -1.0000 -0.0000 -0.0000 0.0000 0 -0.9714 -0.6361 -0.0000 0.0000 0 1.0000 -0.5028 0.5028 0 1.0000 -0.0395 0 0.5156 0.3377 -1.0000 -0.0000 1.0000
za = 6×6
Inf 0 0 0 0 0 0 -0.0172 0 0 0 0 0 0 -0.0580 0 0 0 0 0 0 Inf 0 0 0 0 0 0 Inf 0 0 0 0 0 0 Inf
ma = 6×6
0.0000 0.3750 0.3750 0 0 0 0.0000 -0.4847 0.7401 -0.0000 0 0 0 0.6361 -0.9714 0 0 0 0.0000 -1.0000 -1.0000 0 0 0 1.0000 -0.0009 -0.0316 0 1.0000 0 0.0409 0.0059 -0.0695 -1.0000 0 1.0000
Here we can see the generalized eigenvectors associated to the two zeros and for the equilibrium point .
The columns 2,3 of the matrices , are respectively the two right and left generalized eigenvectos; while the first 4 elements of the 2° and 3° columns of are the zeros state direction and the last two elements of the 2° and 3° columns of and are respectively the zero input and output directions , .
The following are the same quantities for the equilibrium point
[nb,zb,mb] = eig(G_nmp,M)
nb = 6×6
0 0 0 0 0.0000 0 0 0 0 0 0 0.0000 0 1.0000 -1.0000 0.0000 0 -0.0000 0 -0.9716 -0.7740 0 0.0000 0.0000 1.0000 -0.5316 0.5316 -0.0557 1.0000 0.0946 0 0.4953 0.3946 -1.0000 -0.0000 1.0000
zb = 6×6
Inf 0 0 0 0 0 0 0.0128 0 0 0 0 0 0 -0.0562 0 0 0 0 0 0 Inf 0 0 0 0 0 0 Inf 0 0 0 0 0 0 Inf
mb = 6×6
0.0000 0.6755 0.5381 0 0 0 0.0000 -1.0000 1.0000 -0.0000 0 0 0.0000 0.4508 -0.4508 0 0 0 0 -0.5824 -0.4639 0 0 0 1.0000 0.0386 -0.0435 0 1.0000 0 0.0481 -0.0474 -0.0906 -1.0000 0 1.0000
Let us now extract the velues specified before.
csi_z_mp = na(1:4, 2:3) % zero state directions for the two MP zeros of P-
csi_z_mp = 4×2
0 0 0 0 1.0000 -1.0000 -0.9714 -0.6361
v_z_mp = na(5:6, 2:3) % zero input directions for the two MP zeros of P-
v_z_mp = 2×2
-0.5028 0.5028 0.5156 0.3377
csi_z_nmp = nb(1:4, 2:3) % zero input state directions for the two zeros (one MP one NMP) of P+
csi_z_nmp = 4×2
0 0 0 0 1.0000 -1.0000 -0.9716 -0.7740
v_z_nmp = nb(5:6, 2:3) % zero input directions for the two zeros (one MP one NMP) of P+
v_z_nmp = 2×2
-0.5316 0.5316 0.4953 0.3946
w_z_mp = ma(5:6, 2:3) % zero output directions for the two MP zeros of P-
w_z_mp = 2×2
-0.0009 -0.0316 0.0059 -0.0695
w_z_nmp = mb(5:6, 2:3) % zero output directions for the two zeros (one MP one NMP) of P+
w_z_nmp = 2×2
0.0386 -0.0435 -0.0474 -0.0906
eta_z_nmp = mb(1:4, 2:3) % zero output state directions for the two zeros (one MP one NMP) of P+
eta_z_nmp = 4×2
0.6755 0.5381 -1.0000 1.0000 0.4508 -0.4508 -0.5824 -0.4639
From the inspection of the component of , we can see that the direction of the NMP zeros is not a unit vector, so the effect of the RHP zero is distributed between both outputs.
Now we want to show the blocking property of these two pairs of transmition zeros.
In other words, we want to show that the outputs of the two system exept for some small errors due to numerical digit approximations (always present in Matlab), provided that the inputs are set as and and the inizial states .
This means that every couple of inputs with such directions , with the form of an exponential with grow rate equal to , will not be reproduced in non of the outputs provided that the system starts to evolve from the zero state direction , thanks to the presence of the associated transmition zeros .
From the inspection of the component of , we can state both for the MP and the NMP case that, in order to the exhibition by the system of the blocking property of the zeros, the distances of the levels of the first two tanks and from the height at the two equilibrium point and have to be equal to zero, the remaining 3 and 4 have to differ by the quantity specified in the other 2 elements of , and the difference between the current voltages provided and the ones at the two equilibria has to be a signal with exponential form, grow rate equal to , and the proportions in which the exponential should be present at the corresponding input equal to the components of .
Let us try to visalize all of these things by simulations.
t = 1:0.001:1000; % inizializing the time span
 
u1_mp = [zeros(length(t),1), zeros(length(t),1)]; % inizializing u_1 at P_
u2_mp = [zeros(length(t),1), zeros(length(t),1)]; % inizializing u_2 at P_
 
u1_nmp = [zeros(length(t),1), zeros(length(t),1)]; % inizializing u_1 at P+
u2_nmp = [zeros(length(t),1), zeros(length(t),1)]; % inizializing u_2 at P+
 
% building the exponential signals with directions v_z :
 
for i =1 : length(t)
u1_mp(i,1) = exp(zi_mp(2)*t(i))* v_z_mp(1,1);
u1_mp(i,2) = exp(zi_mp(2)*t(i))* v_z_mp(2,1);
end
 
for i =1 : length(t)
u2_mp(i,1) = exp(zi_mp(1)*t(i))* v_z_mp(1,2);
u2_mp(i,2) = exp(zi_mp(1)*t(i))* v_z_mp(2,2);
end
 
for i =1 : length(t)
u1_nmp(i,1) = exp(zi_nmp(2)*t(i))* v_z_nmp(1,1);
u1_nmp(i,2) = exp(zi_nmp(2)*t(i))* v_z_nmp(2,1);
end
 
for i =1 : length(t)
u2_nmp(i,1) = exp(zi_nmp(1)*t(i))* v_z_nmp(1,2);
u2_nmp(i,2) = exp(zi_nmp(1)*t(i))* v_z_nmp(2,2);
end
 
 
% simulating from inizial conditions csi_z:
 
y1 = lsim(G_sys_mp, u1_mp,t, csi_z_mp(:,1));
y2 = lsim(G_sys_mp, u2_mp,t, csi_z_mp(:,2));
 
y3 = lsim(G_sys_nmp, u1_nmp,t, csi_z_nmp(:,1));
y4 = lsim(G_sys_nmp, u2_nmp,t, csi_z_nmp(:,2));
 
figure()
 
subplot(211)
plot(t, y1, 'LineWidth',2); grid
title("MP case: blocking property of the zero z1 = -0.0172");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
 
subplot(212)
plot(t, y2, 'LineWidth',2); grid
title("MP case: blocking property of the zero z2 = -0.0580");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','East')
figure()
subplot(211)
plot(t, y3, 'LineWidth',2); grid
title("NMP case : blocking property of the zero z1 = -0.0562");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
 
subplot(212)
plot(t, y4, 'LineWidth',2); grid
title("NMP case : blocking property of the zero z2 = 0.0128");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','NorthEast')
As we expected, the outputs are identically = 0 except for some small errors of the order of due to the numerical approximations.
In order investigate more in which way the effect of the RPH zero is distributed between the two outputs, we can also compute symbolically the zero output directions through the transfer function, as done in the Johansson's paper.
We start putting :
syms wz1 wz2
 
wz = [wz1; wz2];
 
res = wz'*G
res = 
Then we solve for as suggested in the Johansson's paper.
eq7 = res(2) == 0
eq7 = 
solution = solve(eq7,g2) % solving for gamma2
solution = 
Then we substitute this values in the original expression, we solve for we simplify for .
res1= eval(subs(res, g2, solution)) % substituting gamma2 in the original expression
res1 = 
res2 = collect(res1, wz1) % regrouping terms
res2 = 
eq8 = res2(1) ==0;
s2 = solve(eq8, wz1) % solving for w_z1
s2 = 
eq9 = wz1/wz2 == s2/wz2 % obtaining the Johansson's relation
eq9 = 
This relation shows that as decreases, increases and the effect of the zero is more relevant in the first output. Viceversa, as increases decreases and the zero has a more relevant effect for This is intuitive: as much we provide less flow to one tank, harder to control is its level (for the increased effect of the NMP zero).
Let us now find the symbolical pole directions, remembering that are the input poles directions and the output poles directions, where , are respectively the left and the right eigenvectors associated with the eigenvalues of the matrix A.
[U, Di] = eig(A) % computation of right eigenvectors
U = 
Di = 
[V, Di] = eig(A.') % computation of left eigenvectors (computed as columns)
V = 
Di = 
input_pole_vec = (V.')*B % each row is the input direction associate to i-th eigenvalue (i-th pole)
input_pole_vec = 
output_pole_vec = C*U % each column is the output direction associate to i-th eigenvalue (i-th pole)
output_pole_vec = 
It is interesting to evaluate these quantities in the parameters combination we analyzed in the point 1, the one in which We saw that in such a case we have two zero-pole coincidences in
We know that such a coincidence will bring to a zero-pole cancellation conditionally to several condition.
In particular we have:
So let us evaluate the system numerically with and compute the poles and the zeros directions..
A_canc= eval(subs(A, [A1, A2, A3,A4, t1, t2, t3 ,t4], [A1_n, A2_n, A3_n, A4_n, t1_mp, t2_mp, t3_mp, t4_mp]));
B_canc = eval(subs(B, [A1, A2, A3,A4, g1, g2, k1 , k2], [A1_n, A2_n, A3_n, A4_n, 1, 1, k1_mp, k2_mp]));
C_canc = eval(subs(C, kc, kc_n));
D_canc = eye(2)*0;
 
[U_canc, D] = eig(A_canc) % computation of right eigenvectors
U_canc = 4×4
1.0000 0 -0.8503 0 0 1.0000 0 -0.8315 0 0 0.5263 0 0 0 0 0.5555
D = 4×4
-0.0159 0 0 0 0 -0.0111 0 0 0 0 -0.0419 0 0 0 0 -0.0333
[V_canc, D] = eig(A_canc.') % computation of left eigenvectors
V_canc = 4×4
0 0.5263 0 0 0 0 0 0.5555 1.0000 0.8503 0 0 0 0 1.0000 0.8315
D = 4×4
-0.0419 0 0 0 0 -0.0159 0 0 0 0 -0.0333 0 0 0 0 -0.0111
input_pole_vec_canc = (V_canc.')*B_canc % computation of input poles directions
input_pole_vec_canc = 4×2
0 0 0.0626 0 0 0 0 0.0581
output_pole_vec_canc = C_canc*U_canc % computation of output poles directions
output_pole_vec_canc = 2×4
0.5000 0 -0.4251 0 0 0.5000 0 -0.4158
Gs_canc = [(k1_mp*n1_mp/(1+t1_mp*s)), 0; 0, (k2_mp*n2_mp)/(1+t2_mp*s)] % transfer function
Gs_canc = 
G_ss_canc = ss(A_canc, B_canc, C_canc, D_canc);
tzero(G_ss_canc)
ans = 2×1
-0.0419 -0.0333
As we expected we have two zeros/pole coincidences in .
G_canc = [A_canc,B_canc; C_canc, D_canc]; % inizialization of G (gamma1 = gamma2 = 1)
 
[n_canc,z_canc,m_canc] = eig(G_canc,M) % computing n and m
n_canc = 6×6
0 0 0 0 0.0000 0 0 0 0 0 0 0.0000 0 0 1.0000 0 0 0 0 0 0 1.0000 0 0 1.0000 0 -0.3520 0 1.0000 0 0 1.0000 0 -0.3185 0 1.0000
z_canc = 6×6
Inf 0 0 0 0 0 0 Inf 0 0 0 0 0 0 -0.0419 0 0 0 0 0 0 -0.0333 0 0 0 0 0 0 Inf 0 0 0 0 0 0 Inf
m_canc = 6×6
0.0000 0 0 0 0 0 0 0.0000 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 1.0000 0 0 1.0000 0 0 0 1.0000 0 0 1.0000 0 0 0 1.0000
v_z_canc = n_canc(5:6, 3:4) % zero input directions (gamma1 = gamma2 = 1)
v_z_canc = 2×2
-0.3520 0 0 -0.3185
csi_z_nmp = n_canc(1:4, 3:4) % zero state input directions (gamma1 = gamma2 = 1)
csi_z_nmp = 4×2
0 0 0 0 1 0 0 1
eta_z_canc = m_canc(1:4, 3:4) % zero state output directions (gamma1 = gamma2 = 1)
eta_z_canc = 4×2
0 0 0 0 1 0 0 1
w_z_canc = m_canc(5:6, 3:4) % zero output directions (gamma1 = gamma2 = 1)
w_z_canc = 2×2
0 0 0 0
As we can see, , so with dont's have a pole/zero cancellation with loss of observability. Additionally, , so we have not a zero/pole cancellation with loss of controllability. So this proves that altought there are two zero/pole coincidences as , there are never cancellations and hidden dynamics in the process. This result is coherent with the fact that the realization {A,B,C,D} used is minimal.

3. Limitations due to the presence of NMP zeros.

The first limitation we want to explore is the fact that if G(s) has a NMP-zero at z with output direction , then for internal stability of the closed-loop system it must be guaranteed the validity of the two following interpolation constraints :
i. e. T must have a NMP-zero in the same direction as G(s), and S(z) has an eigenvalue of 1 corresponding to the left eigenvector .
Let's compute a simple stabilizing controller using norm approach, withouth taking into account other things like specifications on the control effort or decoupling aims.
A diagonal weight seems to be sufficient.
s = tf('s');
A1 = 1E-4;
M1 = 1.5;
omb1 = 2.5;
 
wp1 = (s/M1+ omb1)/ (s+omb1*A1);
wp = (s/M1+ omb1)/ (s+omb1*A1)* eye(2); % diagonal weight
 
[KHinf1, Ghinf1, gopt1] = mixsyn(Gs_nmp, wp, [], []);
 
Cloop1 = loopsens(Gs_nmp, KHinf1);
 
r_unit = [ones(length(t),1),ones(length(t),1)]; % dummy reference to check if the controller works
 
figure()
 
ycl1 = lsim(Cloop1.To, r_unit, t);
plot(t, ycl1, 'LineWidth', 2), grid
title("System's steady state with simple H_{inf} controller");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','NorthEast')
xlim([0 200]);
The closed loop system is stable. We can see that the interpolations conditions hold.
%T = (Gs_nmp)/(eye(2)+Gs_nmp); % inizializing complementary sensitivity matrix
T = Cloop1.To;
S = Cloop1.So; % inizializing sensitivity matrix
 
T_z = evalfr(T, zi_nmp(2)); % computing at z = 0.128
S_z = evalfr(S, zi_nmp(2)); % computing at z = 0.128
 
w_z_nmp(:,1)' * T_z % showing the limitation
ans = 1×2
10-12 ×
0.0780 -0.1353
w_z_nmp(:,1)' * S_z % showing the limitation
ans = 1×2
0.0386 -0.0474
We expect then a NMP zero in in T(s) exhibiting as reasonable the blocking property, exactly like G(s).
So let us solve the generalized eigevectors problem in order to compute the zeros input and output directions, and the zero state directions.
T_sys = ss(T);
 
GT = [T_sys.A, T_sys.B; T_sys.C, T_sys.D]; % G matrix for T(s)
MT = [eye(14), T_sys.B*0; T_sys.C*0, T_sys.D*0]; % M matrix for T(s)
 
[nT,zT,mT] = eig(GT,MT)
nT = 16×16
0 0.0000 0.0000 0.0000 0.3665 -0.0798 0.1505 0.0000 0.0000 0.0000 -0.2490 -0.0000 -0.0623 -0.0000 -0.0000 -0.0000 0 0 -0.0000 -0.0000 0.9575 -0.0804 0.0000 -0.7002 0.0000 0.7002 0.0000 -0.0000 0.0000 -0.0009 0.0000 0.0000 0 0 0 -0.0000 -0.2661 -0.0983 -0.0000 1.0000 -0.0000 -1.0000 -0.0000 0.0000 -0.0000 0.0013 0.0000 0.0000 0 0 0.0000 0.0000 1.0000 0.0495 0.1155 -0.0000 -0.0000 0.0000 -0.1911 -0.0000 -0.0478 -0.0000 -0.0000 -0.0000 0 0 0 0.0000 -0.2779 0.0605 -0.1141 -0.0000 -0.0000 -0.0000 0.1887 0.0000 0.0472 0.0000 -0.0000 -0.0000 0 0 0.0000 0.0000 0.2428 0.0897 0.0000 -0.9125 0.0000 0.9125 0.0000 -0.0000 0.0000 -0.0012 -0.0000 -0.0000 0 -0.0000 -0.0000 -0.0000 0.0008 -0.7790 1.0000 0.0000 1.0000 0.3740 1.0000 0.6955 1.0000 -0.7131 0.0000 0.0000 0 -0.0000 -0.0000 0.0000 0.0003 1.0000 -0.9740 -0.0000 -0.8936 -0.5244 -0.9740 -1.0000 -0.9740 1.0000 -0.0000 -0.0000 0 -0.0000 -0.0000 -0.0000 -0.3665 0.0798 -0.0025 -0.0000 0.0000 -0.0000 -0.1820 0.0000 0.0549 -0.0000 0.0000 0.0000 0 0.0000 0.0000 0.0000 -0.9575 0.0804 -0.0000 0.7002 -0.0000 0.5802 -0.0000 0.1943 -0.0000 -0.0328 -0.0000 -0.0000
zT = 16×16
Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Inf 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0562 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0128 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0158 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0109 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0256 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.0109 0 0 0 0 0 0
mT = 16×16
0.0000 -0.0000 -0.0000 -0.0000 -0.7423 0.5369 1.0000 0.3650 0 0 0 0 0 0 0 0 0.0000 -0.0000 -0.0000 -0.0000 0.3712 -0.2685 -0.4619 0.0633 -0.0000 0.1232 -0.5000 1.0000 0.4513 0.8767 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.6526 -0.7138 -0.3826 0.0722 -0.0000 0.1406 -0.4141 0.7002 0.3738 1.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.5000 0.2881 -0.0964 0.5121 -0.9876 -0.0000 0.0482 0.0000 0.5779 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.4730 1.0000 -0.1582 1.0000 -1.0000 -0.0000 0.0791 0.0000 0.9481 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -1.0000 -0.5762 -0.0648 -0.5314 0.0000 1.0000 -0.0702 0.0000 0.0633 0.4399 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0762 0.0000 -0.0000 0.0000 -1.0000 0.0000 0.9026 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 0.2458 -0.0000 0.1232 0.0000 1.0000 -0.0000 0.8767 0.0000 0.0000
zT(6,6) % the NMP zeros present also in G(s)
ans = 0.0128
csi_zT = nT(1:14, 6); % zero input state directions
v_zT = nT (15:16, 6); % zero input directions
We can verify the presence of the zero we expected.
All this means that for every controller such that it is able to stabilize the plant G(s) and the resulting closed loop system exhibits the stability property, we will have a blocking NMP zero in 0.0128 In T(s).
Remembering that, in the case of no presence of disturbances or noises,
we can state that there will be a lack of control autority with respect to the references of the form and the inizial states , sinche the outputs with respect to such references will be identically equal to 0. Let's visualize it with a simulation.
 
rT = [ones(length(t),1),ones(length(t),1)]; % inizializing the 2 references
 
for i =1 : length(t)
rT(i,1) = exp(zi_nmp(2)*t(i))* v_zT(1,1); % imposing values and directions
rT(i,2) = exp(zi_nmp(2)*t(i))* v_zT(2,1);
end
 
y5 = lsim(T_sys,rT,t, csi_zT(1:14)); % outputs
 
figure()
plot(t, y5, 'LineWidth',2); grid
title("Limitation on T: blocking property of the zero z_{2} = 0.0128");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','NorthEast')
We can see that in such directions the reference can not be reproduced at steady state: the 2 outputs are always identically equal to zero (except for the problems of approximation that brings initially the outputs to reach values of the order of ).
An additional limitation is derived from the fact that in the SISO design we require:
In a multivariable context, this is equivalent to require:
where is the largest singular value of S(s).
SInce , we obtain an imposition of a bandwith limitation equal to : in the case of z real, like this case .
So, if we want to require a sensitivity peak to be not more than 2 as usual, and a very small error at steady state (A << 1), we obtain that the bandiwidth can not be greater than .
So in our design we have to take into account the fact that
We can see from similations that only requiring a bandiwidth below this quantity the MATLAB function mixsyn() is able to find a stabilizing controller coherent with our specifications (.
If we try to require more bandiwith than this, we obtain a controller such that again the closed loop bandiwidth for S(s) is not greater than , and consequently an higher γ, increasing as it increases the quantity of bandiwith required. (Remembering that γ is an index of how much the resulting controller is far from the point of view of the performances from our initial specifications).
A2 = 1E-4;
M2 = 2;
omb2 = 0.0064;
 
om = logspace(-2,5,1000);
 
wp2 = (s/M2+ omb2)/ (s+omb2*A2);
wp22 = (s/M2+ omb2)/ (s+omb2*A2)* eye(2); % diagonal weight
 
[KHinf2, Ghinf2, gamma2 gopt2] = mixsyn(Gs_nmp, wp22, [], []);
 
gamma2 % we have a gamma very close to 1
gamma2 = 1.0040
Cloop2 = loopsens(Gs_nmp, KHinf2);
 
r_unit2 = [ones(length(t),1),ones(length(t),1)]; % dummy reference to check if the controller works
 
figure()
 
ycl2 = lsim(Cloop2.To, r_unit2, t);
plot(t, ycl2, 'LineWidth', 2), grid
title("System's steady state with simple H_{inf} controller (NMP)");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
xlim([0 1000]);
ylim([-0.5 1.3])
sigmaWS2 = sigma(wp22,om);
sigmaWS2Inv = sigma(1/wp22,om);
sigmaSo2 = sigma(Cloop2.So,om);
 
figure()
semilogx(om,1./sigmaWS2(2,:),om, sigmaWS2Inv(1,:),'--',... % plotting the minimum singular value of the weight , the max of the inverted weight, and the max of S(s)
om, sigmaSo2(1,:),'LineWidth',2), grid;
legend('1/s_{min}(wS2)','s_{max}(wS2^{-1})','s_{max}(S_o)','Location','northeast');
xlim([1E-2,1E+2]);
ylim([0 3])
ylabel('dB','FontName','courier','FontSize',12)
xlabel('$\omega$','Interpreter','latex','FontSize',14)
We can see that requiring exactly = 0.0064 we obtain a controller perfectly coherent with our specifications, such that the maximum singular value of the closed loop sensitivity function S(s) is perfectly upper bounded from the minimum singular value of the designed weight (or of the maximum singular value of the inverse of our weight).

4. Robustness Analysis

Minimum Phase case

Let us now consider some robustness issues, taking as a case of study a situation in which the parameters are uncertain and vary in intervals between +10% and -10% of their nominal value. Initially we start to consider the system G(s) exhibiting Minimum Phase Properties.
NMatlabRed = [0.8500 0.3250 0.0980];
NMatlabYellow = [0.929 0.694 0.125 ];
NMatlabBlue = [0 0.4470 0.7410];
NMatlabViolet = [0.4940 0.1840 0.5560];
NMatlabGreen = [0.4660 0.6740 0.1880];
NMatlabCyan = [0.3010 0.7450 0.9330];
NMatlabBordeaux = [0.6350 0.0780 0.1840];
 
% uncertain k_i
 
k1_mp_u = ureal('K1',3.33,'Range',[(3.33-(3.33/10)), (3.33+(3.33/10))]); % cm^3/Vs
k2_mp_u = ureal('K2',3.35,'Range',[(3.35-(3.35/10)), (3.35+(3.35/10))]); % cm^3/Vs
 
% uncertain gamma_i
 
g1_mp_u = ureal('G1',0.7,'Range',[(0.7-(0.7/10)), (0.7+(0.7/10))]);
g2_mp_u = ureal('G2',0.6,'Range',[(0.6-(0.6/10)), (0.6+(0.6/10))]);
 
Gs_mp_unc = [(g1_mp_u*k1_mp_u*n1_mp/(1+t1_mp*s)), ((1-g2_mp_u)*k2_mp_u*n1_mp)/ ((1+t3_mp*s)*(1+t1_mp*s)); ((1-g1_mp_u)*k1_mp_u*n2_mp)/ ((1+t2_mp*s)*(1+t4_mp*s)), (g2_mp_u*k2_mp_u*n2_mp)/(1+t2_mp*s)]
Uncertain continuous-time state-space model with 2 outputs, 2 inputs, 6 states. The model uncertainty consists of the following blocks: G1: Uncertain real, nominal = 0.7, range = [0.63,0.77], 2 occurrences G2: Uncertain real, nominal = 0.6, range = [0.54,0.66], 2 occurrences K1: Uncertain real, nominal = 3.33, range = [3,3.66], 2 occurrences K2: Uncertain real, nominal = 3.35, range = [3.02,3.69], 2 occurrences Model Properties Type "Gs_mp_unc.NominalValue" to see the nominal value and "Gs_mp_unc.Uncertainty" to interact with the uncertain elements.
We want to represent this parametric uncertainty as an unstructured input multiplicative uncertainty of the form:
with and thw weight such that it satisfy the requirement that:
, with the set of all possible perturbed plants.
We compute a 3° order diagonal weight through the command ucover().
w =logspace (-3,2, 1000); % frequency span
 
figure()
 
Gs_mp_unc_sampled = gridureal(Gs_mp_unc,30);
G_star_mp = inv(Gs_mp)* (Gs_mp_unc_sampled- Gs_mp); % showing the singular values of the matrix G^-1( G_p - G)
sigma(G_star_mp, w);
hold on;
 
[PwM1,InfoM1] = ucover(Gs_mp_unc_sampled,Gs_mp,3, 'InputMult'); % computing a suitable weight
 
[m, p] = bode(InfoM1.W1, w);
plot_weight = semilogx(w, 20*log10(m(1,:)));
set(plot_weight,'LineWidth',4, 'Color',NMatlabRed), grid on;
Leg2 = legend('$\sigma$($G^{-1}$($G_{u}$-G)','$w_m$ with n = 3','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
grid on;
title('Weight for singular values, multiplicative uncertanity','FontSize',10);
Let us visualize the resulting input multiplicative representation of the uncertainty in terms of relative error:
figure()
 
[PwM1,InfoM1] = ucover(Gs_mp_unc_sampled,Gs_mp,[3 3], 'InputMult'); % choiche of a diagonal weight
Rel_P = (Gs_mp- Gs_mp_unc_sampled )/Gs_mp; % relative error
bodemag(Rel_P,'b--',InfoM1.W1,'r'), grid on;
Leg2 = legend('sampled','$w_m$ with n = 3','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
grid on;
title('Relative error, multiplicative uncertanity','FontSize',10);
figure
bodemag(PwM1,'b',Gs_mp,'r'), grid on;
Leg2 = legend('uncertain','nominal','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
grid on
title('multiplicative uncertanity representation','FontSize',10);
We can see from this plots that the most uncertain situation is at low frequency, while at high frequency the diagonal elements of the uncertainty representations are both very close to the nominal process behavior .This result is also highlighted in the Vadigepalli, Gatzke, Doyle's paper.
Let us now compute two different stabilizing decentralized controllers for the nominal case. We want to compare the performance of two different controllers in terms of reference reproduction's capabilities and also in terms of robust stability properties.
The first controller is the decentralized controller PI proposed in the Johansson's paper, the second one is another diagonal controller designed with an norm control approach.
% controller for the G(s) with MP property
 
K1_mp = 3.0;
T1_mp = 30;
 
K2_mp = 2.7;
T2_mp = 40;
 
PI1_mp = K1_mp*(1+ (1/(T1_mp*s)));
PI2_mp = K2_mp*(1+ (1/(T2_mp*s)));
 
PI_mp = [ PI1_mp ,0; 0, PI2_mp] % decentralized PI controller with Johansson's parameters
PI_mp = From input 1 to output... 90 s + 3 1: -------- 30 s 2: 0 From input 2 to output... 1: 0 108 s + 2.7 2: ----------- 40 s Continuous-time transfer function. Model Properties
Cloop3 = loopsens(Gs_mp, PI_mp);
 
r_unit2 = [ones(length(t),1),ones(length(t),1)];
 
figure()
 
ycl3 = lsim(Cloop3.To, r_unit2, t);
plot(t, ycl3, 'LineWidth', 2), grid
title("System's steady state with Decentralized PI controller (MP)");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
xlim([0 300]);
ylim([-0.5 1.3])
On the other hand, we need for comparison an controller for the MP case. We will use the input-output paring derived to be better in the Relative Array Gain analysis:
. This implies that the performance weight of the decentralized controller will have off-diagonal components and the diagonal ones = 0.
Since we don't have costraints on the control effort, we just design a weight for the sensitivity function S(s).
A3 = 1E-3;
M3 = 1.2;
omb3 = 0.2;
 
om = logspace(-2,5,1000);
 
wp3 = (s/M3+ omb3)/ (s+omb3*A3);
wp33 = (s/M3+ omb3)/ (s+omb3*A3)* [0,1;1,0]; % off-diagonal weight
 
[KHinf3, Ghinf3, gamma3 gopt2] = mixsyn(Gs_mp, wp33, [], []);
 
gamma3 % we have a gamma very close to 1
gamma3 = 0.8363
Cloop5 = loopsens(Gs_mp, KHinf3);
 
figure()
 
ycl5 = lsim(Cloop5.To, r_unit2, t);
plot(t, ycl5, 'LineWidth', 2), grid
title("System's steady state with H_{inf} controller (MP)");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
xlim([0 30]);
ylim([-0.5 1.3])
We can see the differences in terms of nominal performance: the controller allows a faster transient and not overshoot.
Let us check if this 2 controllers are robust stable, i.e. if they satisfy the necessary and sufficient condition .
L_MP_PI = PI_mp*PwM1; % open loop function PI
L_MP_Hinf = KHinf3*PwM1; % open loop function H-infinity
 
CL_MP_PI = feedback(L_MP_PI, eye(2)); % closed loop function PI
CL_MP_Hinf= feedback(L_MP_Hinf, eye(2)); % closed loop function H-infinity
 
[stabmarg1,wcu1] = robstab(CL_MP_PI)
stabmarg1 = struct with fields:
LowerBound: 1.7088 UpperBound: 1.7123 CriticalFrequency: 1.0000e-04
wcu1 = struct with fields:
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
[stabmarg2,wcu2] = robstab(CL_MP_Hinf)
stabmarg2 = struct with fields:
LowerBound: 1.7102 UpperBound: 1.7137 CriticalFrequency: 1.0000e-04
wcu2 = struct with fields:
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
Being the upper and lower bounds striclty >1 in both case, we can state that both controllers are able to robust stabilize the uncertain system with MP property.
Let us check if they can also guarantee robust performance, i.d. if where M is the component of the augmented plant the only possible source of instability of the overall control structured system and is the structured singular value.
[perfmarg_MP_Hinf,wcu,report,info1] = robustperf(CL_MP_Hinf)
perfmarg_MP_Hinf = struct with fields:
LowerBound: 0.9642 UpperBound: 0.9643 CriticalFrequency: 0.0310
wcu = struct with fields:
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
report = 5×129 char array
'Uncertain system does not achieve performance robustness to modeled uncertainty. ' ' -- The tradeoff of model uncertainty and system gain is balanced at a level of 96.4% of the modeled uncertainty. ' ' -- A model uncertainty of 96.4% can lead to input/output gain of 1.04 at 0.031 rad/seconds. ' ' -- Sensitivity with respect to each uncertain element is: ' ' 13% for Gs_mp_unc_sampled_InputMultDelta. Increasing Gs_mp_unc_sampled_InputMultDelta by 25% decreases the margin by 3.25%.'
info1 = struct with fields:
Sensitivity: [1×1 struct] Frequency: [113×1 double] BadUncertainValues: [113×1 struct] MussvBnds: [1×2 frd] MussvInfo: [1×1 struct]
[perfmarg_MP_PI,wcu,report,info2] = robustperf(CL_MP_PI)
perfmarg_MP_PI = struct with fields:
LowerBound: 0.8010 UpperBound: 0.8022 CriticalFrequency: 0.0567
wcu = struct with fields:
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
report = 5×129 char array
'Uncertain system does not achieve performance robustness to modeled uncertainty. ' ' -- The tradeoff of model uncertainty and system gain is balanced at a level of 80.1% of the modeled uncertainty. ' ' -- A model uncertainty of 80.2% can lead to input/output gain of 1.25 at 0.0567 rad/seconds. ' ' -- Sensitivity with respect to each uncertain element is: ' ' 15% for Gs_mp_unc_sampled_InputMultDelta. Increasing Gs_mp_unc_sampled_InputMultDelta by 25% decreases the margin by 3.75%.'
info2 = struct with fields:
Sensitivity: [1×1 struct] Frequency: [69×1 double] BadUncertainValues: [69×1 struct] MussvBnds: [1×2 frd] MussvInfo: [1×1 struct]
While we had robust stability for both control systems, we don't have robust performance capabilities, being the lower and the upper bounds for robust performance strictly <1 in both cases.
Let us visualize how the structured singulare value varies with variation of frequencies.
% extracting the values for the H inf controller
 
a = info1.MussvBnds(1,1);
b = info1.MussvBnds(1,2);
max_mp_Hinf = zeros(length(info1.Frequency),1);
for i=1:size(info1.Frequency)
max_mp_Hinf(i)= a.Response(:,:,i);
end
min_mp_Hinf = zeros(length(info1.Frequency),1);
for i=1:length(info1.Frequency)
min_mp_Hinf(i)= b.Response(:,:,i);
end
 
% extracting the values for the PI controller
 
c = info2.MussvBnds(1,1);
d = info2.MussvBnds(1,2);
max_mp_PI = zeros(length(info2.Frequency),1);
for i=1:length(info2.Frequency)
max_mp_PI(i)= c.Response(:,:,i);
end
min_mp_PI = zeros(length(info2.Frequency),1);
for i=1:length(info2.Frequency)
min_mp_PI(i)= d.Response(:,:,i);
end
 
figure()
 
subplot(211)
semilogx(info1.Frequency, max_mp_Hinf,'b--',info1.Frequency, min_mp_Hinf,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds H-inf CL');
xlim([1.e-4,1.e+2])
ylim([0,1.5])
grid on;
 
subplot(212)
semilogx(info2.Frequency, max_mp_PI,'b--',info2.Frequency, min_mp_PI,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds PI CL');
ylim([0,1.5])
xlim([1.e-4,1.e+2])
grid on;
This plot confirms our expectations: both controllers can not achieve robust performance since the structured singular value is >1 for some ω, but we can see that the one is significantly closer to achieve RP: the upper and lower bounds are higher than the ones of the PI closed loop and also the condition < 1 is vialoted at very few frequencies.
Let us try to design a controller able to guarantee robust performance trhough μ-synthesis control.
We have to build the augumented plant and run the D-K algorithm.
opts = hinfsynOptions;
opts.Display = 'on';
 
% weights for S(s)
 
A4 = 1E-2;
M4 = 1.28;
omb4 = 40;
 
wp44 = ((s/M4+ omb4)/ (s+omb4*A4))* [0,1;1,0]; % off diagonal weight
 
% inizializing the inputs and the outputs of each block
 
PwM1.u = 'u'; % control input
PwM1.y = 'y'; % not perturbed output
wp44.u = 'v';
wp44.y = 'z'; % performance variable
 
Sum3 = sumblk('v = r - y', 2); % measured variables
 
Pe_MP = connect(PwM1,wp44,Sum3, {'r','u'},{'z','v' }); % augmented plant
 
ncont = 2; % two control inputs
nmeas = 2; % two measurements
[Krob1,rpMU1] = musyn(Pe_MP ,nmeas,ncont)
D-K ITERATION SUMMARY: ----------------------------------------------------------------- Robust performance Fit order ----------------------------------------------------------------- Iter K Step Peak MU D Fit D 1 30.04 30.04 30.28 8 2 10.88 10.88 10.88 8 3 3.278 2.964 2.983 20 4 2.784 1.797 1.807 12 5 1.704 1.454 1.46 12 6 1.422 1.301 1.306 12 7 1.28 1.204 1.212 20 8 1.2 1.144 1.153 20 9 1.135 1.094 1.103 20 10 1.08 1.05 1.059 20 Best achieved robust performance: 1.05 Krob1 = A = x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32 x33 x34 x1 -0.01328 1 6.232e-10 2.358e-10 1.329e-10 3.849e-14 -1.685e-11 -6.507e-10 -2.576e-10 -1.384e-10 -1.029e-09 2.088e-10 1.229e-09 -2.482e-10 -9.192e-10 1.019e-09 -1.712e-10 -6.92e-11 -2.254e-10 1.551e-10 6.101e-11 2.05e-10 1.15e-08 -1.446e-08 -1.529e-14 6.513e-12 2.574e-10 9.181e-11 5.624e-11 1.411e-14 -6.451e-12 -2.522e-10 -9.273e-11 -5.493e-11 x2 -0.000435 0.1163 7.053 3.388 1.283 0.0003546 -0.08741 -3.66 -1.222 -0.8192 -10.87 2.009 6.158 -2.254 -9.522 4.464 -1.551 -0.6273 -2.051 1.507 0.5947 1.99 196.5 -0.09352 -0.0001713 0.06759 2.63 1.096 0.5263 0.0001599 -0.04579 -1.889 -0.6306 -0.4227 x3 -0.01262 4.378 171.7 114.5 43.34 0.01198 -2.953 -123.6 -41.29 -27.68 -367.3 67.87 208 -76.16 -321.7 150.8 -52.41 -21.19 -69.28 50.92 20.09 67.24 6640 -3.159 -0.005786 2.284 88.86 37.03 17.78 0.005402 -1.547 -63.8 -21.31 -14.28 x4 -3.859e-18 9.36e-16 3.563e-14 -104.5 86.75 4.058e-18 1.101e-15 3.026e-14 2.207e-14 4.621e-15 -1.871e-13 2.507e-14 -9.62e-14 -2.316e-14 -1.559e-13 -1.075e-13 -1.54e-14 -6.227e-15 -2.059e-14 1.977e-14 7.848e-15 2.604e-14 3.636e-12 4.29e-12 -2.484e-18 4.505e-16 1.888e-14 8.151e-15 3.787e-15 2.709e-18 -1.044e-16 -8.243e-15 -2.272e-16 -2.299e-15 x5 -0.01034 3.589 140.8 -37.4 -64.38 0.00982 -2.42 -101.3 -33.84 -22.68 -301.1 55.63 170.5 -62.43 -263.7 123.6 -42.96 -17.37 -56.79 41.74 16.47 55.11 5442 -2.59 -0.004742 1.872 72.83 30.36 14.57 0.004428 -1.268 -52.3 -17.46 -11.71 x6 2.58e-13 -1.132e-10 -4.484e-09 -1.532e-09 -9.953e-10 -0.01328 1 6.224e-09 2.694e-09 1.282e-09 4.474e-09 -1.203e-09 -1.265e-08 1.581e-09 4.252e-09 -1.112e-08 1.103e-09 4.46e-10 1.445e-09 -8.668e-10 -3.394e-10 -1.148e-09 -1.468e-08 2.209e-07 7.324e-14 -4.507e-11 -1.757e-09 -5.942e-10 -3.916e-10 -5.808e-14 5.246e-11 1.956e-09 7.832e-10 4.145e-10 x7 0.0002309 -0.09807 -3.906 -1.304 -0.8743 -0.0002789 0.1444 7.787 3.75 1.494 3.458 -1.1 -11.67 1.429 3.392 -10.27 0.9931 0.4015 1.301 -0.7984 -0.3126 -1.057 -0.08686 241.9 6.595e-05 -0.03976 -1.554 -0.519 -0.3479 -5.961e-05 0.06113 2.235 0.989 0.4477 x8 0.007799 -3.313 -132 -44.06 -29.54 -0.007343 5.326 196.5 126.7 50.47 116.8 -37.15 -394.2 48.29 114.6 -346.9 33.55 13.56 43.96 -26.97 -10.56 -35.72 -2.934 8172 0.002228 -1.343 -52.51 -17.54 -11.75 -0.002014 2.065 75.51 33.41 15.12 x9 5.032e-18 -2.445e-15 -9.731e-14 -3.06e-14 -2.223e-14 -4.398e-18 4.421e-15 1.601e-13 -104.5 86.75 4.948e-14 -2.123e-14 -3.363e-13 3.098e-14 5.391e-14 -3.029e-13 2.186e-14 8.836e-15 2.851e-14 -1.476e-14 -5.747e-15 -1.96e-14 7.871e-13 6.695e-12 9.919e-19 -9.492e-16 -3.659e-14 -1.182e-14 -8.286e-15 -5.543e-19 1.237e-15 4.475e-14 1.888e-14 9.308e-15 x10 0.006393 -2.716 -108.2 -36.12 -24.21 -0.006019 4.366 161.1 -27.37 -58.54 95.77 -30.45 -323.1 39.58 93.94 -284.3 27.5 11.12 36.03 -22.11 -8.657 -29.28 -2.405 6699 0.001826 -1.101 -43.04 -14.37 -9.634 -0.001651 1.693 61.89 27.39 12.4 x11 -0.02447 8.644 338.4 118.9 96.35 0.02375 -5.985 -250 -83.47 -55.95 -756.6 133.6 412.3 -148.6 -658.5 299.5 -101.6 -41.07 -134.2 101 39.82 133.3 1.373e+04 -6.46 -0.01077 3.657 146.2 50.06 37.07 0.01058 -3.098 -127.5 -42.58 -28.55 x12 -0.02446 8.644 338.4 118.9 96.35 0.02375 -5.985 -250 -83.47 -55.95 -756.6 133.5 412.3 -148.6 -658.4 299.5 -101.6 -41.06 -134.2 100.9 39.82 133.3 1.373e+04 -6.678 -0.01077 3.657 146.2 50.06 37.07 0.01058 -3.098 -127.5 -42.58 -28.54 x13 -8.713e-07 0.0002955 0.01154 0.004915 0.002331 8.55e-07 -0.0001538 -0.00688 -0.00195 -0.001603 -0.02889 0.0205 0.01 -0.005283 -0.02494 0.006176 -0.003599 -0.001455 -0.004764 0.003707 0.001464 0.004894 0.4412 0.1244 -4.065e-07 0.0001266 0.005086 0.001922 0.001085 4.071e-07 -9.919e-05 -0.004197 -0.001328 -0.0009527 x14 0.01472 -6.638 -263.3 -87.92 -58.94 -0.01324 10.48 384.8 142.8 109.7 225.9 -66.9 -791.2 90.22 220 -701.3 63.19 25.54 82.68 -47.79 -18.69 -63.34 -5.998 1.68e+04 0.00384 -2.625 -102 -34.06 -22.83 -0.002866 3.157 116.7 41.98 31.08 x15 9.9e-07 -0.0004812 -0.01915 -0.00602 -0.004376 -8.652e-07 0.0008701 0.03152 0.0142 0.006394 0.009721 -0.004175 -0.0662 0.03735 0.0106 -0.05963 0.004301 0.001738 0.005609 -0.002904 -0.00113 -0.003855 0.1553 1.318 1.95e-07 -0.0001868 -0.0072 -0.002325 -0.00163 -1.088e-07 0.0002434 0.008806 0.003716 0.001832 x16 0.01472 -6.639 -263.3 -87.93 -58.94 -0.01324 10.48 384.9 142.8 109.7 225.9 -66.91 -791.3 90.29 220 -701.5 63.19 25.55 82.69 -47.8 -18.69 -63.34 -5.783 1.681e+04 0.003841 -2.625 -102 -34.07 -22.83 -0.002867 3.157 116.7 41.99 31.09 x17 0.000617 0.1222 3.429 1.765 0.9 0.0005514 -0.4215 -16.32 -5.451 -3.652 -79.14 0.9272 9.275 1.744 -60.04 8.126 2.664 1.081 3.619 2.346 0.9055 3.121 1543 -1.011 0.001259 -1.8 -61.53 -49.43 4.797 -3.235e-05 -0.1402 -5.188 -1.733 -1.161 x18 -0.0003107 -0.06033 -1.685 -0.87 -0.4439 -0.0002735 0.2103 8.14 2.72 1.822 39.49 -0.4511 -4.611 -0.8856 29.96 -4.046 -1.328 -0.5698 -1.822 -1.163 -0.449 -1.548 -770.2 0.8034 -0.0006303 0.9004 30.78 24.72 -2.396 1.722e-05 0.06986 2.583 0.8633 0.5782 x19 7.805e-05 0.01496 0.4162 0.2154 0.11 6.802e-05 -0.0525 -2.031 -0.6789 -0.4546 -9.862 0.1108 1.147 0.2236 -7.48 1.007 0.3336 0.1494 0.4321 0.2893 0.1117 0.385 192.4 -0.2097 0.0001577 -0.2251 -7.696 -6.181 0.5985 -4.465e-06 -0.01742 -0.644 -0.2152 -0.1441 x20 -0.001277 -0.3061 -9.79 -3.271 -2.192 0.002551 0.07231 -1.456 1.614 -0.4688 -7.173 13.09 -39.49 -9.731 -10.49 -45.98 -5.745 -2.325 -7.762 11.15 4.421 14.72 -0.3622 1697 -0.00116 0.02189 2.183 0.7279 0.4886 0.002374 -1.959 -68.9 -51.61 3.107 x21 0.0006396 0.1525 4.872 1.629 1.091 -0.001276 -0.03503 0.7682 -0.7886 0.2425 3.592 -6.549 19.66 4.873 5.253 22.91 2.877 1.164 3.887 -5.566 -2.237 -7.363 0.4968 -846.8 0.0005804 -0.01116 -1.1 -0.3666 -0.2462 -0.001187 0.9799 34.46 25.81 -1.551 x22 -0.0001599 -0.03809 -1.217 -0.4069 -0.2723 0.0003191 0.008674 -0.195 0.1957 -0.06121 -0.8968 1.637 -4.908 -1.218 -1.312 -5.722 -0.7195 -0.2912 -0.9721 1.392 0.5655 1.816 -0.1778 211.5 -0.0001451 0.0028 0.2754 0.09175 0.06164 0.0002968 -0.245 -8.615 -6.452 0.3876 x23 -8.047e-16 3.318e-13 1.31e-11 4.716e-12 2.852e-12 7.467e-16 -4.237e-13 -1.595e-11 -6.654e-12 -3.332e-12 -1.732e-11 3.95e-12 3.144e-11 -4.917e-12 -1.585e-11 2.701e-11 -3.41e-12 -1.379e-12 -4.479e-12 2.894e-12 1.136e-12 3.829e-12 -0.4 -4.761e-10 -2.674e-16 1.343e-13 5.273e-12 1.833e-12 1.163e-12 2.327e-16 -1.445e-13 -5.512e-12 -2.121e-12 -1.183e-12 x24 2.499e-16 -1.029e-13 -4.063e-12 -1.464e-12 -8.84e-13 -2.319e-16 1.309e-13 4.932e-12 2.055e-12 1.03e-12 5.4e-12 -1.228e-12 -9.714e-12 1.527e-12 4.937e-12 -8.339e-12 1.059e-12 4.279e-13 1.391e-12 -8.999e-13 -3.533e-13 -1.191e-12 -4.455e-11 -0.4 8.326e-17 -4.166e-14 -1.636e-12 -5.689e-13 -3.608e-13 -7.258e-17 4.476e-14 1.708e-12 6.565e-13 3.667e-13 x25 -0.00146 0.4749 18.5 8.174 3.665 0.001447 -0.1614 -8.245 -1.608 -2.052 -51.59 8.352 9.167 -8.838 -44.23 3.064 -6 -2.426 -7.953 6.392 2.526 8.435 829.9 410.4 -0.014 1.206 8.323 3.201 1.763 0.0007309 -0.1481 -6.473 -1.92 -1.492 x26 -0.0001797 0.07268 2.796 1.188 0.5746 0.0002135 -0.0533 -2.256 -0.6986 -0.5149 -9.072 1.14 2.932 -1.156 -7.575 2.003 -0.7393 -0.2986 -0.9762 0.9192 0.3622 1.214 147.6 19.61 -0.0001121 -0.0428 1.134 0.4325 0.2431 9.017e-05 -0.02814 -1.157 -0.3748 -0.2611 x27 0.0002828 0.2757 9.298 4.213 2.127 0.0009954 -0.5629 -22.09 -7.38 -4.944 -99.66 3.198 16.97 -0.4936 -76.99 13.86 1.308 0.5409 1.796 4.233 1.648 5.616 1879 -1.572 0.001214 -1.928 -65.6 0.5666 0.2325 0.0001365 -0.2053 -7.805 -2.608 -1.747 x28 -0.0009465 0.3345 13.1 5.38 2.692 0.0009189 -0.2316 -9.675 -3.231 -2.165 -29.23 5.167 15.96 -5.748 -25.44 11.6 -3.929 -1.589 -5.193 3.906 1.541 5.158 418.3 -0.4068 -0.0004168 0.1415 5.656 -102.4 87.96 0.0004093 -0.1199 -4.934 -1.648 -1.105 x29 0.003771 -1.025 -41.35 -16.66 -8.324 -0.00262 0.4047 18.07 6.032 4.045 27.63 -16.7 -45.76 21.09 32.04 -31.99 15.77 6.384 20.89 -11.14 -4.41 -14.69 -24.05 -0.04842 0.002554 -2.109 -74.88 -138.6 -104.3 -0.001419 0.28 12.05 4.024 2.698 x30 0.001913 -0.8858 -35.17 -11.5 -7.932 -0.001704 1.471 53.72 23.8 10.97 25.8 -8.483 -111.3 11.75 25.86 -99.26 8.248 3.334 10.78 -6.01 -2.346 -7.968 100.7 2104 0.0004579 -0.3479 -13.48 -4.451 -3.03 -0.0136 1.431 15.79 6.512 3.312 x31 5.275e-05 -0.05776 -2.224 -0.6867 -0.5111 2.08e-06 0.09425 3.262 1.56 0.6531 0.3984 0.04637 -8.304 0.265 0.442 -7.861 0.2255 0.09104 0.2842 0.1043 0.04374 0.1342 23.18 199.8 -2.269e-05 -0.01765 -0.6337 -0.1999 -0.1447 1.007e-05 -0.05235 0.6119 0.3038 0.1207 x32 -0.00119 -0.4359 -14.7 -4.912 -3.29 0.002617 0.2347 4.087 4.276 0.6501 -4.548 13.44 -55.24 -9.395 -8.297 -61.07 -5.4 -2.186 -7.336 11.61 4.614 15.28 -0.9056 2085 -0.001223 -0.01475 0.8966 0.2979 0.2006 0.002576 -2.114 -74.33 -1.965 -1.785 x33 0.0004843 -0.2185 -8.665 -2.894 -1.94 -0.0004356 0.3448 12.66 5.553 2.596 7.44 -2.201 -26.01 2.971 7.245 -23.05 2.079 0.8405 2.721 -1.573 -0.6148 -2.084 -0.359 475.1 0.0001264 -0.08637 -3.357 -1.121 -0.7513 -9.435e-05 0.1039 3.838 -102.9 87.56 x34 -0.002861 0.4932 21.69 7.239 4.854 0.00384 -1.15 -45.96 -18.12 -9.577 -32.66 19.58 56 -19.26 -34.98 39.72 -12.52 -5.063 -16.6 15.64 6.174 20.64 0.01414 -141.6 -0.001494 0.3241 13.8 4.608 3.089 0.002478 -2.137 -75.82 -138.9 -104.5 B = u1 u2 x1 -9.5e-10 9.238e-10 x2 -4.467e-11 2.183e-11 x3 -4.602e-12 -1.894e-12 x4 -3.006e-13 -2.744e-13 x5 -1.354e-13 4.559e-14 x6 1.205e-09 -1.413e-08 x7 -2.223e-11 -2.847e-11 x8 -1.218e-13 1.665e-12 x9 -6.529e-14 -4.28e-13 x10 -8.181e-13 -5.564e-13 x11 -99.36 0.002303 x12 -99.29 0.01625 x13 -0.03647 -0.007969 x14 0.005347 -88.57 x15 -0.01288 -0.08427 x16 -0.01248 -88.71 x17 -217.4 0.03289 x18 108.6 -0.03552 x19 -27.14 0.009449 x20 -0.00229 -194.2 x21 -0.02502 96.96 x22 0.01069 -24.23 x23 8 3.044e-11 x24 3.677e-12 8 x25 -68.61 -26.27 x26 -15.14 -1.258 x27 -254.5 0.05914 x28 -34.57 0.01106 x29 -79.35 0.02508 x30 -8.399 -134.5 x31 -1.922 -15.57 x32 0.03168 -227.8 x33 0.01269 -30.38 x34 0.02958 -68.38 C = x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32 x33 x34 y1 -0.09898 34.35 1347 472.3 383.8 0.09399 -23.17 -970 -323.9 -217.1 -2882 532.5 1632 -597.6 -2524 1183 -411.2 -166.3 -543.6 399.5 157.6 527.5 5.209e+04 -24.79 -0.04539 17.92 697.2 290.6 139.5 0.04238 -12.14 -500.6 -167.2 -112.1 y2 0.06119 -25.99 -1035 -345.7 -231.8 -0.05761 41.79 1542 568.4 439.7 916.7 -291.5 -3093 378.9 899.2 -2722 263.2 106.4 344.9 -211.6 -82.86 -280.3 -23.02 6.412e+04 0.01748 -10.54 -412 -137.6 -92.22 -0.0158 16.2 592.4 262.1 118.7 D = u1 u2 y1 0 0 y2 0 0 Continuous-time state-space model. Model Properties
rpMU1 = 1.0502
After 10 iterations, the best attempt of musyn() using D-K algorithm ends with a order controller able to guarantee a maximum structured singular value of 1.052, so not robust performant in absolute terms, but very close.
We expect lower and upper bounds for robust performance indices very close to one.
L_MP_rob = Krob1*PwM1; % open loop function robust performance controller (MP)
 
CL_MP_rob = feedback(L_MP_rob, eye(2)); % closed loop robust performance controller (MP)
 
[perfmarg_MP_rob,wcu,report,info3] = robustperf(CL_MP_rob)
perfmarg_MP_rob = struct with fields:
LowerBound: 0.9805 UpperBound: 0.9806 CriticalFrequency: 12.3801
wcu = struct with fields:
Gs_mp_unc_sampled_InputMultDelta: [2×2 ss]
report = 5×128 char array
'Uncertain system does not achieve performance robustness to modeled uncertainty. ' ' -- The tradeoff of model uncertainty and system gain is balanced at a level of 98.1% of the modeled uncertainty. ' ' -- A model uncertainty of 98.1% can lead to input/output gain of 1.02 at 12.4 rad/seconds. ' ' -- Sensitivity with respect to each uncertain element is: ' ' 14% for Gs_mp_unc_sampled_InputMultDelta. Increasing Gs_mp_unc_sampled_InputMultDelta by 25% decreases the margin by 3.5%.'
info3 = struct with fields:
Sensitivity: [1×1 struct] Frequency: [121×1 double] BadUncertainValues: [121×1 struct] MussvBnds: [1×2 frd] MussvInfo: [1×1 struct]
Let us visualize the behavior of the structured singular value in a suitable range of frequencies.
% extracting the structured singular values for the µ-syn robust controller
 
p = info3.MussvBnds(1,1);
o = info3.MussvBnds(1,2);
max_mp_rob = zeros(length(info3.Frequency),1);
 
for i=1:length(info3.Frequency)
max_mp_rob(i)= p.Response(:,:,i);
end
 
min_mp_rob = zeros(length(info3.Frequency),1);
 
for i=1:length(info3.Frequency)
min_mp_rob(i)= o.Response(:,:,i);
end
 
figure()
 
semilogx(info3.Frequency, max_mp_rob,'b--',info3.Frequency, min_mp_rob,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds robust performant CL');
xlim([1.e-4,1.e+6])
ylim([0,1.5])
grid on;
We can see that µ is always <1 exept for a small neighborood of the critical frequency . This result is reasonably the best achieved in terms of robust performance,
better than the nominal controller which achieves a robust performance index around 96.42%, and than the PI that achieves about 80%.
Let us compute the step response and try to visualize some possible behaviors (20) of closed loop system under the µ-syn controller.
L_MP_rob = Krob1* PwM1; % robust open loop
CL_MP_rob = feedback(L_MP_rob, eye(2)); % robust closed loop
 
figure()
samples_MP = usample(CL_MP_rob ,20);
lsim(samples_MP , 'b',r_unit2, t);
title("System's steady state with µ-synthesis controller (MP)");
xlabel("t");
ylabel("y(t)");
xlim([0 20]);
grid on;

Non-Minimum Phase case

% uncertain k_i
 
k1_nmp_u = ureal('K1',3.14,'Range',[(3.14-(3.14/10)), (3.14+(3.14/10))]); % cm^3/Vs
k2_nmp_u = ureal('K2',3.29,'Range',[(3.29-(3.29/10)), (3.29+(3.29/10))]); % cm^3/Vs
 
% uncertain gamma_i
 
g1_nmp_u = ureal('G1',0.43,'Range',[(0.43-(0.43/10)), (0.43+(0.43/10))]);
g2_nmp_u = ureal('G2',0.34,'Range',[(0.34-(0.34/10)), (0.34+(0.34/10))]);
 
Gs_nmp_unc = [(g1_nmp_u*k1_nmp_u*n1_nmp/(1+t1_nmp*s)), ((1-g2_nmp_u)*k2_nmp_u*n1_nmp)/ ((1+t3_nmp*s)*(1+t1_nmp*s)); ((1-g1_nmp_u)*k1_nmp_u*n2_nmp)/ ((1+t2_nmp*s)*(1+t4_nmp*s)), (g2_nmp_u*k2_nmp_u*n2_nmp)/(1+t2_nmp*s)]
Uncertain continuous-time state-space model with 2 outputs, 2 inputs, 6 states. The model uncertainty consists of the following blocks: G1: Uncertain real, nominal = 0.43, range = [0.387,0.473], 2 occurrences G2: Uncertain real, nominal = 0.34, range = [0.306,0.374], 2 occurrences K1: Uncertain real, nominal = 3.14, range = [2.83,3.45], 2 occurrences K2: Uncertain real, nominal = 3.29, range = [2.96,3.62], 2 occurrences Model Properties Type "Gs_nmp_unc.NominalValue" to see the nominal value and "Gs_nmp_unc.Uncertainty" to interact with the uncertain elements.
We want to represent this parametric uncertainty as an unstructured input multiplicative uncertainty of the form:
with and thw weight such that it satisfy the requirement that:
, with the set of all possible perturbed plants.
We compute a 3° order diagonal weight through the command ucover().
w =logspace (-3,2, 1000); % frequency span
 
figure()
 
Gs_nmp_unc_sampled = gridureal(Gs_nmp_unc,30);
G_star_nmp = inv(Gs_nmp)* (Gs_nmp_unc_sampled- Gs_nmp); % showing the singular values of the matrix G^-1( G_p - G)
sigma(G_star_nmp, w);
hold on;
 
[PwM2,InfoM2] = ucover(Gs_nmp_unc_sampled,Gs_nmp,3, 'InputMult'); % computing a suitable weight
 
[m, p] = bode(InfoM2.W1, w);
plot_weight = semilogx(w, 20*log10(m(1,:)));
set(plot_weight,'LineWidth',4, 'Color',NMatlabRed), grid on;
Leg2 = legend('$\sigma$($G^{-1}$($G_{u}$-G)','$w_m$ with n = 3','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
grid on;
title('Weight for singular values, multiplicative uncertanity','FontSize',10);
Let us visualize the resulting input multiplicative representation of the uncertainty in terms of relative error:
figure()
 
[PwM2,InfoM2] = ucover(Gs_nmp_unc_sampled,Gs_nmp,[3 3], 'InputMult'); % choiche of a diagonal weight
Rel_P = (Gs_nmp- Gs_nmp_unc_sampled )/Gs_nmp; % relative error
bodemag(Rel_P,'b--',InfoM2.W1,'r'), grid on;
Leg2 = legend('sampled','$w_m$ with n = 3','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
grid on;
title('Relative error, multiplicative uncertanity','FontSize',10);
figure
bodemag(PwM2,'b',Gs_nmp,'r'), grid on;
Leg2 = legend('uncertain','nominal','Interpreter','latex');
set(Leg2,'FontSize',12,'Location','SouthEast','Interpreter','latex');
grid on
title('multiplicative uncertanity representation','FontSize',10);
We can see from this plots that the most uncertain situation is at low frequency, while at high frequency the diagonal elements of the uncertainty representations are both very close to the nominal process behavior .This result is also highlighted in the Vadigepalli, Gatzke, Doyle's paper.
Let us now compute two different stabilizing decentralized controllers for the nominal case.
The first controller is the decentralized controller PI proposed in the Johansson's paper, the second one is another diagonal controller designed with an norm control approach.
% controller for the G(s) with NMP property
 
K1_nmp = 1.5;
T1_nmp = 110;
 
K2_nmp = -0.12;
T2_nmp = 220;
 
PI1_nmp = K1_nmp*(1+ (1/(T1_nmp*s)));
PI2_nmp = K2_nmp*(1+ (1/(T2_nmp*s)));
 
PI_nmp = [PI1_nmp, 0; 0,PI2_nmp] % decentralized PI controller with Johansson's parameters
PI_nmp = From input 1 to output... 165 s + 1.5 1: ----------- 110 s 2: 0 From input 2 to output... 1: 0 -26.4 s - 0.12 2: -------------- 220 s Continuous-time transfer function. Model Properties
 
Cloop4 = loopsens(Gs_nmp, PI_nmp);
t_nmp = 1:0.001:4000; % inizializing the time span for the non minimum phase transient!
r_unit3 = [ones(length(t_nmp),1),ones(length(t_nmp),1)];
 
figure()
 
ycl4 = lsim(Cloop4.To, r_unit3, t_nmp);
plot(t_nmp, ycl4, 'LineWidth', 2), grid
title("System's steady state with Decentralized PI controller (NMP)");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
xlim([0 4000]);
ylim([-0.5 1.9])
Note the difference in terms of time scale: the PI controller produces an output of the system with NMP property with a transient about 10 times slower than in the MP case.
For the controller we can use the one designed in the previous section for the NMP case, the one that takes into account the bandwith limitations imposed by the presence of the transmittion zero.
%weight selection
A7 = 1E-4;
M7 = 2;
omb7 = 0.0064;
 
 
om = logspace(-2,5,1000);
 
wp55 = (s/M7+ omb7)/ (s+omb7*A7)* [0,1;1,0]; % off-diagonal weight
 
[KHinf7, Ghinf7, gamma7, gopt7] = mixsyn(Gs_nmp, wp55, [], []);
 
gamma7 % we have a gamma very close to 1
gamma7 = 1.0040
Cloop7 = loopsens(Gs_nmp, KHinf7);
 
figure()
 
ycl7 = lsim(Cloop7.To, r_unit2, t);
plot(t, ycl7, 'LineWidth', 2), grid
title("System's steady state with H_{inf} controller (MP)");
xlabel("t");
ylabel("y(t)");
legend('y1(t)','y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
xlim([0 500]);
ylim([-0.5 1.3])
We can see the differences in terms of nominal performance: the controller allows a faster transient and not overshoot.
In both situations the outputs are slower than the corresponding MP case of about a order of magnitude.
Let us check if this 2 controllers are robust stable, i.e. if they satisfy the necessary and sufficient condition .
L_NMP_PI = PI_nmp*PwM2; % open loop function PI
L_NMP_Hinf = KHinf7*PwM2; % open loop function H-infinity
 
CL_NMP_PI = feedback(L_NMP_PI, eye(2)); % closed loop function PI
CL_NMP_Hinf= feedback(L_NMP_Hinf, eye(2)); % closed loop function H-infinity
 
[stabmarg1,wcu1] = robstab(CL_NMP_PI)
stabmarg1 = struct with fields:
LowerBound: 1.7651 UpperBound: 1.7687 CriticalFrequency: 0.0052
wcu1 = struct with fields:
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
[stabmarg2,wcu2] = robstab(CL_NMP_Hinf)
stabmarg2 = struct with fields:
LowerBound: 3.6340 UpperBound: 3.6413 CriticalFrequency: 8.9079e-04
wcu2 = struct with fields:
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
Being the upper and lower bounds striclty >1 in both case, we can state that both controllers are able to robust stabilize the uncertain system with MP property.
Let us check if they can also guarantee robust performance, i.d. if where M is the component of the augmented plant the only possible source of instability of the overall control structured system and is the structured singular value.
[perfmarg_MP_Hinf,wcu,report,info1] = robustperf(CL_NMP_Hinf)
perfmarg_MP_Hinf = struct with fields:
LowerBound: 0.6966 UpperBound: 0.6966 CriticalFrequency: 0.0647
wcu = struct with fields:
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
report = 5×131 char array
'Uncertain system does not achieve performance robustness to modeled uncertainty. ' ' -- The tradeoff of model uncertainty and system gain is balanced at a level of 69.7% of the modeled uncertainty. ' ' -- A model uncertainty of 69.7% can lead to input/output gain of 1.44 at 0.0647 rad/seconds. ' ' -- Sensitivity with respect to each uncertain element is: ' ' 23% for Gs_nmp_unc_sampled_InputMultDelta. Increasing Gs_nmp_unc_sampled_InputMultDelta by 25% decreases the margin by 5.75%.'
info1 = struct with fields:
Sensitivity: [1×1 struct] Frequency: [100×1 double] BadUncertainValues: [100×1 struct] MussvBnds: [1×2 frd] MussvInfo: [1×1 struct]
[perfmarg_MP_PI,wcu,report,info2] = robustperf(CL_NMP_PI)
perfmarg_MP_PI = struct with fields:
LowerBound: 0.3708 UpperBound: 0.3709 CriticalFrequency: 0.0066
wcu = struct with fields:
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
report = 5×128 char array
'Uncertain system does not achieve performance robustness to modeled uncertainty. ' ' -- The tradeoff of model uncertainty and system gain is balanced at a level of 37.1% of the modeled uncertainty. ' ' -- A model uncertainty of 37.1% can lead to input/output gain of 2.7 at 0.00658 rad/seconds. ' ' -- Sensitivity with respect to each uncertain element is: ' ' 20% for Gs_nmp_unc_sampled_InputMultDelta. Increasing Gs_nmp_unc_sampled_InputMultDelta by 25% decreases the margin by 5%.'
info2 = struct with fields:
Sensitivity: [1×1 struct] Frequency: [87×1 double] BadUncertainValues: [87×1 struct] MussvBnds: [1×2 frd] MussvInfo: [1×1 struct]
As in the MP case, while we had robust stability for both control systems, we don't have robust performance capabilities, being the lower and the upper bounds for robust performance strictly <1 in both cases.
Let us visualize how the structured singulare value varies with variation of frequencies.
% extracting the values for the H inf controller
 
a = info1.MussvBnds(1,1);
b = info1.MussvBnds(1,2);
max_nmp_Hinf = zeros(length(info1.Frequency),1);
for i=1:size(info1.Frequency)
max_nmp_Hinf(i)= a.Response(:,:,i);
end
min_nmp_Hinf = zeros(length(info1.Frequency),1);
for i=1:length(info1.Frequency)
min_nmp_Hinf(i)= b.Response(:,:,i);
end
 
% extracting the values for the PI controller
 
c = info2.MussvBnds(1,1);
d = info2.MussvBnds(1,2);
max_nmp_PI = zeros(length(info2.Frequency),1);
for i=1:length(info2.Frequency)
max_nmp_PI(i)= c.Response(:,:,i);
end
min_nmp_PI = zeros(length(info2.Frequency),1);
for i=1:length(info2.Frequency)
min_nmp_PI(i)= d.Response(:,:,i);
end
 
figure()
 
subplot(211)
semilogx(info1.Frequency, max_nmp_Hinf,'b--',info1.Frequency, min_nmp_Hinf,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds H-inf CL');
xlim([1.e-4,1.e+2])
ylim([0,2.5])
grid on;
 
subplot(212)
semilogx(info2.Frequency, max_nmp_PI,'b--',info2.Frequency, min_nmp_PI,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds PI CL');
ylim([0,3])
xlim([1.e-4,1.e+2])
grid on;
Here we have another important limitation of the NMP case: the robust performance condition is much more violated by both controllers than the MP case. Anyaway the again looks like to behave better than the PI.
Let us try to design a controller able to guarantee robust performance trhough μ-synthesis control.
We have to build the augumented plant and run the D-K algorithm: the requirements on the weights have to be coherent with the limitations imposed by the NMP property of the system.
opts = hinfsynOptions;
opts.Display = 'on';
 
% weights for S(s)
 
A8 = 1E-1;
M8 = 2;
omb8 = 0.0064; % maximal bandiwith allowed for NMP zero
 
wp66 = ((s/M8+ omb8)/ (s+omb8*A8))* [0,1;1,0]; % off diagonal weight
 
% inizializing the inputs and the outputs of each block
 
PwM2.u = 'u'; % control input
PwM2.y = 'y'; % not perturbed output
wp66.u = 'v';
wp66.y = 'z'; % performance variable
 
Sum4 = sumblk('v = r - y', 2); % measured variables
 
Pe_NMP = connect(PwM2,wp66,Sum4, {'r','u'},{'z','v' }); % augmented plant
 
ncont = 2; % two control inputs
nmeas = 2; % two measurements
[Krob2,rpMU2] = musyn(Pe_NMP ,nmeas,ncont)
D-K ITERATION SUMMARY: ----------------------------------------------------------------- Robust performance Fit order ----------------------------------------------------------------- Iter K Step Peak MU D Fit D 1 1.522 1.391 1.405 4 2 1.265 1.222 1.234 12 3 1.206 1.19 1.202 12 4 1.199 1.186 1.198 12 5 1.2 1.188 1.2 12 Best achieved robust performance: 1.19 Krob2 = A = x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x1 -0.01107 0.23 0.1125 -0.006708 -0.08067 -0.04269 0.4043 -0.3527 -0.8182 0.3497 0.8875 -0.6259 -0.01516 -0.01827 0.01828 0.01098 0.01482 -0.0108 -0.4932 0.7911 0.004226 0.05411 0.02667 -0.001738 -0.02093 -0.01107 x2 -1.89e-16 -0.03478 1 1.753e-16 2.156e-15 1.137e-15 -1.462e-14 6.226e-15 1.774e-14 -6.438e-15 -2.581e-14 1.466e-14 2.434e-16 2.944e-16 -2.94e-16 -2.4e-16 -3.243e-16 2.339e-16 1.273e-14 -1.464e-14 -4.721e-17 -5.763e-16 -3.043e-16 4.399e-17 5.403e-16 2.849e-16 x3 0.07354 0.6774 -1.074 -0.09991 -1.202 -0.6359 6.022 -5.254 -12.19 5.209 13.22 -9.323 -0.2258 -0.2721 0.2723 0.1635 0.2208 -0.1608 -7.346 11.78 0.06295 0.8059 0.3972 -0.02588 -0.3117 -0.1649 x4 -0.006885 -0.07939 -0.04219 -0.01148 0.2242 0.1124 -0.6562 0.3644 0.741 -0.3956 -1.266 0.5324 0.01256 0.01527 -0.01519 -0.01657 -0.02239 0.01598 0.7574 -0.6578 -0.001834 -0.02137 -0.01134 0.004165 0.05318 0.02682 x5 1.029e-16 1.255e-15 6.625e-16 -9.538e-17 -0.03478 1 7.943e-15 -3.388e-15 -9.661e-15 3.504e-15 1.403e-14 -7.987e-15 -1.325e-16 -1.602e-16 1.6e-16 1.306e-16 1.764e-16 -1.273e-16 -6.925e-15 7.969e-15 2.57e-17 3.137e-16 1.657e-16 -2.393e-17 -2.939e-16 -1.55e-16 x6 -0.1026 -1.182 -0.6284 0.06744 0.591 -1.075 -9.773 5.428 11.04 -5.893 -18.86 7.93 0.1871 0.2274 -0.2263 -0.2469 -0.3335 0.2381 11.28 -9.799 -0.02732 -0.3182 -0.1689 0.06203 0.7922 0.3996 x7 0.04902 0.4652 0.5305 -0.06672 -0.8049 -0.4257 4.136 -3.353 -8.043 3.34 8.858 -6.242 -0.1401 -0.1693 0.1692 0.1068 0.1442 -0.1049 -4.895 7.583 0.01464 0.1512 0.1398 -0.01721 -0.2078 -0.1099 x8 0.0756 0.6569 0.9166 -0.1127 -1.354 -0.7167 6.576 -5.996 -13.99 5.915 14.66 -10.74 -0.2512 -0.3037 0.3036 0.1851 0.2499 -0.182 -8.282 13.43 0.02367 0.234 0.2435 -0.02921 -0.3516 -0.186 x9 0.02788 0.3401 0.1796 -0.02583 -0.3176 -0.1674 2.147 -0.9025 -2.62 0.9493 3.794 -2.167 -0.0359 -0.04342 0.04336 0.03538 0.0478 -0.03448 -1.876 2.16 0.006966 0.08503 0.0449 -0.006481 -0.07961 -0.04198 x10 -0.114 -1.312 -0.6976 0.06679 0.5354 0.9045 -11.1 6.103 12.26 -6.671 -21.38 8.75 0.2099 0.2552 -0.2539 -0.2728 -0.3689 0.2638 12.72 -10.95 -0.03043 -0.354 -0.188 0.02215 0.2124 0.2431 x11 -0.01345 -0.164 -0.0866 0.01249 0.1536 0.08097 -1.044 0.4432 1.261 -0.4427 -1.842 1.042 0.01733 0.02096 -0.02093 -0.0171 -0.0231 0.01666 0.9069 -1.042 -0.00336 -0.04101 -0.02166 0.003134 0.03849 0.0203 x12 -0.07143 -0.832 -0.4416 0.04678 0.4323 0.539 -6.667 3.527 7.481 -3.806 -12.65 5.481 0.1235 0.15 -0.1494 -0.1547 -0.2092 0.1498 7.331 -6.589 -0.01882 -0.221 -0.1172 0.01443 0.1475 0.1433 x13 -0.06701 -0.838 -0.4262 0.04465 0.5466 0.2884 -3.979 1.747 4.154 -1.777 -7.049 3.212 -0.05918 0.2008 -0.1379 -0.06365 -0.08597 0.06217 3.263 -3.886 -0.1765 -2.485 -0.8002 0.01128 0.138 0.07282 x14 0.05024 0.6232 0.3209 -0.0378 -0.4635 -0.2445 3.269 -1.424 -3.653 1.458 5.789 -2.911 0.01952 -0.1494 0.09497 0.05302 0.07162 -0.05175 -2.755 3.24 0.09243 1.294 0.4271 -0.009522 -0.1167 -0.06156 x15 -0.01066 -0.1326 -0.06798 0.007676 0.09406 0.04962 -0.6681 0.2932 0.7355 -0.2995 -1.185 0.581 -0.007333 0.0401 -0.0437 -0.01083 -0.01463 0.01057 0.56 -0.6627 -0.02263 -0.3176 -0.1037 0.001936 0.0237 0.01251 x16 0.04815 0.5784 0.3059 -0.06938 -0.8721 -0.4477 3.399 -1.867 -4.919 1.975 6.484 -3.958 -0.0693 -0.08397 0.08374 -0.07588 0.2361 -0.1409 -3.742 4.022 0.01226 0.1476 0.07805 -0.2025 -2.856 -0.9161 x17 -0.03762 -0.4544 -0.2402 0.0473 0.5911 0.3056 -2.758 1.38 3.727 -1.45 -5.107 3.025 0.05211 0.06311 -0.06296 0.03215 -0.1643 0.08726 2.787 -3.06 -0.009513 -0.1151 -0.06085 0.1044 1.467 0.4785 x18 0.00845 0.1019 0.05389 -0.01095 -0.137 -0.07071 0.6177 -0.3137 -0.8414 0.33 1.149 -0.6813 -0.0118 -0.01429 0.01426 -0.009242 0.04517 -0.04356 -0.6329 0.6909 0.00214 0.02587 0.01367 -0.02589 -0.3641 -0.1182 x19 -2.435e-13 -2.97e-12 -1.568e-12 2.258e-13 2.776e-12 1.464e-12 -1.881e-11 8.021e-12 2.287e-11 -8.294e-12 -3.321e-11 1.89e-11 3.136e-13 3.793e-13 -3.788e-13 -3.092e-13 -4.177e-13 3.014e-13 -0.00064 -1.886e-11 -6.084e-14 -7.427e-13 -3.922e-13 5.665e-14 6.959e-13 3.67e-13 x20 6.803e-15 8.298e-14 4.381e-14 -6.308e-15 -7.756e-14 -4.089e-14 5.255e-13 -2.241e-13 -6.388e-13 2.317e-13 9.278e-13 -5.281e-13 -8.761e-15 -1.06e-14 1.058e-14 8.638e-15 1.167e-14 -8.419e-15 -4.58e-13 -0.00064 1.7e-15 2.075e-14 1.096e-14 -1.583e-15 -1.944e-14 -1.025e-14 x21 0.3046 3.714 1.962 -0.2835 -3.485 -1.838 23.52 -10.07 -28.81 10.41 41.57 -23.85 -0.3932 -0.4757 0.4749 0.3881 0.5243 -0.3782 -20.58 23.7 0.04716 0.929 0.4908 -0.07112 -0.8736 -0.4607 x22 0.02298 0.2808 0.1482 -0.02009 -0.2466 -0.13 1.459 -0.741 -2.272 0.7633 2.698 -1.922 -0.02913 -0.03523 0.03518 0.02792 0.03772 -0.02723 -1.46 1.789 0.005727 0.03525 1.037 -0.005053 -0.06197 -0.03269 x23 -0.07915 -0.9912 -0.5037 0.05042 0.6162 0.3252 -4.174 2.027 5.038 -2.057 -7.588 3.965 0.08472 0.1023 -0.1023 -0.0727 -0.09819 0.07106 3.688 -4.574 -0.2122 -2.989 -2.872 0.01276 0.1559 0.0823 x24 -0.2405 -2.934 -1.549 0.2218 2.725 1.438 -18.72 7.91 22.52 -8.177 -32.97 18.61 0.3094 0.3743 -0.3737 -0.3046 -0.4115 0.2969 16.19 -18.61 -0.06009 -0.7336 -0.3873 0.0267 0.6835 0.3607 x25 -0.02375 -0.2893 -0.1528 0.02282 0.2808 0.148 -2.04 0.7928 2.156 -0.8216 -3.523 1.754 0.03088 0.03737 -0.0373 -0.03098 -0.04185 0.03018 1.656 -1.834 -0.005943 -0.07247 -0.03827 0.005717 0.0355 1.037 x26 0.0473 0.5672 0.3 -0.07072 -0.8897 -0.4565 3.62 -1.865 -4.746 1.978 6.777 -3.77 -0.06893 -0.08354 0.0833 0.08063 0.1089 -0.07818 -3.787 3.959 0.01207 0.1451 0.07674 -0.2102 -2.965 -2.861 B = u1 u2 x1 -1.656e-14 2.029e-14 x2 -1.842e-14 2.233e-14 x3 -2.539e-15 3.117e-15 x4 -6.646e-15 6.774e-15 x5 9.966e-15 -1.221e-14 x6 2.487e-16 -3.136e-16 x7 0.5508 -1.76 x8 -1.054 -0.8439 x9 2.68 -3.325 x10 -0.3871 -0.8198 x11 -1.325 1.578 x12 -1.613 1.283 x13 -6.085 0.9564 x14 4.616 -2.502 x15 -0.9605 0.4071 x16 0.8523 -5.454 x17 -1.788 4.296 x18 0.3617 -0.9548 x19 0.0625 2.888e-11 x20 6.597e-13 0.0625 x21 29.19 -36.99 x22 1.016 -3.683 x23 -5.263 2.481 x24 -24.03 28.2 x25 -3.093 2.192 x26 1.767 -4.435 C = x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 y1 0.3471 3.197 3.947 -0.4715 -5.671 -3.001 28.42 -24.79 -57.51 24.58 62.38 -44 -1.065 -1.284 1.285 0.7716 1.042 -0.7589 -34.67 55.61 0.2971 3.803 1.875 -0.1221 -1.471 -0.7784 y2 -0.484 -5.58 -2.966 0.3183 2.789 3.947 -46.12 25.62 52.09 -27.81 -89 37.43 0.8829 1.073 -1.068 -1.165 -1.574 1.123 53.24 -46.24 -0.1289 -1.502 -0.7973 0.2927 3.738 1.886 D = u1 u2 y1 0 0 y2 0 0 Continuous-time state-space model. Model Properties
rpMU2 = 1.1862
After only 5 iterations, the best attempt of musyn() using D-K algorithm ends with a order controller able to guarantee a maximum structured singular value of 1.18, so not robust performant. Also the attempt of relaxing the requirements in the weight selections parameters is not sufficient to achieve better performances: we accepted also 1% of steady state error and an maximal overshoot of 2.
We expect lower and upper bounds for robust performance indices significantly lower than the MP case.
L_NMP_rob = Krob2*PwM2; % open loop function robust performance controller (MP)
 
CL_NMP_rob = feedback(L_NMP_rob, eye(2)); % closed loop robust performance controller (MP)
 
[perfmarg_NMP_rob,wcu,report,info4] = robustperf(CL_NMP_rob)
perfmarg_NMP_rob = struct with fields:
LowerBound: 0.7627 UpperBound: 0.7627 CriticalFrequency: 0.0875
wcu = struct with fields:
Gs_nmp_unc_sampled_InputMultDelta: [2×2 ss]
report = 5×128 char array
'Uncertain system does not achieve performance robustness to modeled uncertainty. ' ' -- The tradeoff of model uncertainty and system gain is balanced at a level of 76.3% of the modeled uncertainty. ' ' -- A model uncertainty of 76.3% can lead to input/output gain of 1.31 at 0.0875 rad/seconds. ' ' -- Sensitivity with respect to each uncertain element is: ' ' 24% for Gs_nmp_unc_sampled_InputMultDelta. Increasing Gs_nmp_unc_sampled_InputMultDelta by 25% decreases the margin by 6%.'
info4 = struct with fields:
Sensitivity: [1×1 struct] Frequency: [89×1 double] BadUncertainValues: [89×1 struct] MussvBnds: [1×2 frd] MussvInfo: [1×1 struct]
Let us visualize the behavior of the structured singular value in a suitable range of frequencies.
% extracting the structured singular values for the µ-syn robust controller
 
p = info4.MussvBnds(1,1);
o = info4.MussvBnds(1,2);
max_nmp_rob = zeros(length(info4.Frequency),1);
 
for i=1:length(info4.Frequency)
max_nmp_rob(i)= p.Response(:,:,i);
end
 
min_nmp_rob = zeros(length(info4.Frequency),1);
 
for i=1:length(info4.Frequency)
min_nmp_rob(i)= o.Response(:,:,i);
end
 
figure()
 
semilogx(info4.Frequency, max_nmp_rob,'b--',info4.Frequency, min_nmp_rob,'r-','LineWidth',1.5);
legend('Lower Bound','Upper Bound');
title('µ bounds robust performant CL');
xlim([1.e-4,1.e+6])
ylim([0,1.5])
grid on;
We can see that µ is not always <1: in a significantly large neighborood of the critical frequency the structured singular value is about 1,27. This result is reasonably the best achieved in terms of robust performance, better than the nominal controller which achieves a robust performance index around 69.66%, and than the PI that achieves about 37,08%, but it is consistently worse than the MP case.
This is another symptom of how much difficult is to control a NMP MIMO system than one will all zeros
Let us compute the step response and try to visualize some possible behaviors (10) of closed loop system under the µ-syn controller.
L_NMP_rob = Krob2* PwM2; % robust open loop
CL_NMP_rob = feedback(L_NMP_rob, eye(2)); % robust closed loop
 
t1= 0:0.01:400;
r_unit4 = [ones(length(t1),1),ones(length(t1),1)];
 
figure()
samples_NMP = usample(CL_NMP_rob ,10);
lsim(samples_NMP , 'b',r_unit4, t1);
title("System's steady state with µ-synthesis controller (NMP)");
xlabel("t");
ylabel("y(t)");
xlim([0 400]);
grid on;
In the end, we can note two intresting aspects:
This concludes our robustness analysis.

5. Decoupling approaches.

Let us try to design two possible decoupler in order to separate the input-output channels minimizing or hopefully nihilating the interactions. We want to find the matrix such that the MIMO system behaves like two decoupled SISO systems, with the possibility of controlling them with two additional SISO controllers (Simplified Decoupling approach), or including the control in the design of (Inverted Decoupling approach). In the first case the overall control system will look like:
apparent_plant.png
while in the second one will substantially coincide with

Simplified Decoupling

The first attempt is based on what it is called "Simplified feedforward decoupling". Matrix is composed in this way:
with and .
and the overall "apparent process" will be:
. Note that the zeros of are the same as .
Furthermore, since both and are non-zero, all transmition zeros not coinciding with poles (i.e. solutions of , will appear in the resulting open loop system as zeros of the elements on the diagonal.
This is an important limitation, since we have a case in which the system's transfer function exhibitis a Non Minimum Phase transmition zero.
So we expect two positive zeros, one on each element of the diagonal, that are more difficult to manage since NMP zeros in the SISO case don't have directions and their effects are present in all situations, while in the MIMO case they act only along the zeros directions.
% minimum phase case
 
d12mp = - (Gs_mp(1,2)/(Gs_mp(1,1))); %element out of the diagonal of D(s)
d21mp= - (Gs_mp(2,1)/(Gs_mp(2,2))); %element out of the diagonal of D(s)
Ds_mp = [1, d12mp; d21mp, 1];
 
L_mp_app= simplify(zpk(Gs_mp*Ds_mp)) % open loop
L_mp_app = From input 1 to output... 0.041625 (s+0.05802) (s+0.01718) (s+0.01595) (s+0.01107) 1: -------------------------------------------------------- (s+0.04186) (s+0.03334) (s+0.01595)^2 (s+0.01107) 6.1876e-20 s (s+0.07544) (s-0.01294) 2: ------------------------------------ (s+0.03334)^2 (s+0.01107)^3 From input 2 to output... -1.585e-19 s (s^2 + 0.001953) 1: ----------------------------- (s+0.04186)^2 (s+0.01595)^3 0.031406 (s+0.05802) (s+0.01718) (s+0.01595) (s+0.01107) 2: -------------------------------------------------------- (s+0.04186) (s+0.03334) (s+0.01595) (s+0.01107)^2 Continuous-time zero/pole/gain model. Model Properties
Since the elements off-diagonal exhibits coefficient of the order of , we can consider them identically equal to zero (they are actual different just for numerical issues for MATLAB). So we can consider our decoupling problem solved.
In order to avoid numerical problems in the next few computations, let's substitute zero off-diagonal instead of the remaining very small values, whose presence is a consequence of numerical MATLAB approximations.
L_mp= [L_mp_app(1,1), 0; 0, L_mp_app(2,2)]; % open looop actually decoupled
 
zero(L_mp(1,1))
ans = 4×1
-0.0580 -0.0172 -0.0159 -0.0111
zero(L_mp(2,2))
ans = 4×1
-0.0580 -0.0172 -0.0159 -0.0111
pole(L_mp(1,1))
ans = 5×1
-0.0419 -0.0333 -0.0159 -0.0159 -0.0111
pole(L_mp(2,2))
ans = 5×1 complex
-0.0419 + 0.0000i -0.0333 + 0.0000i -0.0159 + 0.0000i -0.0111 + 0.0000i -0.0111 - 0.0000i
We can see from the poles and the zeros computation that there are several cancellations in the remaining two SISO elements on the diagonal, but this fact does not constitute a problem since they all are with negative real part. The design of a stabilizing controller can start from here.
Let us explore the Non Minimum Phase Case.
%non minimum phase case
d12nmp = - (Gs_nmp(1,2)/(Gs_nmp(1,1)));
d21nmp= - (Gs_nmp(2,1)/(Gs_nmp(2,2)));
Ds_nmp = zpk([1, d12nmp; d21nmp, 1])
Ds_nmp = From input 1 to output... 1: 1 -0.028515 (s+0.01094) 2: ----------------------- (s+0.01782) (s+0.01094) From input 2 to output... -0.041223 (s+0.01582) 1: ----------------------- (s+0.02563) (s+0.01582) 2: 1 Continuous-time zero/pole/gain model. Model Properties
L_nmp_app = simplify(zpk(Gs_nmp*Ds_nmp))
L_nmp_app = From input 1 to output... 0.024111 (s+0.05623) (s+0.01582) (s+0.01094) (s-0.01278) 1: -------------------------------------------------------- (s+0.01582)^2 (s+0.01782) (s+0.02563) (s+0.01094) 1.4677e-19 (s-0.000503) (s+0.01094) (s+0.02335) 2: ----------------------------------------------- (s+0.01782)^2 (s+0.01094)^3 From input 2 to output... -3.5894e-19 (s+0.01582) (s^2 + 0.0215s + 0.0002711) 1: --------------------------------------------------- (s+0.01582)^3 (s+0.02563)^2 0.017478 (s+0.05623) (s+0.01582) (s+0.01094) (s-0.01278) 2: -------------------------------------------------------- (s+0.01782) (s+0.01582) (s+0.02563) (s+0.01094)^2 Continuous-time zero/pole/gain model. Model Properties
Again we have the same numerical problems as before.
L_nmp= [L_nmp_app(1,1), 0; 0, L_nmp_app(2,2)]; % open looop actually decoupled
zero(L_nmp(1,1))
ans = 4×1
-0.0562 0.0128 -0.0158 -0.0109
zero(L_nmp(2,2))
ans = 4×1
-0.0562 0.0128 -0.0109 -0.0158
pole(L_nmp)
ans = 10×1
-0.0158 -0.0256 -0.0158 -0.0178 -0.0109 -0.0178 -0.0109 -0.0109 -0.0256 -0.0158
We can see the presene of NMP zeros in the two transfer functions, as well as the presence of several cancellations.

Inverted Decoupling

The goal now is to design a controller able to solve the decoupling problem and at the same time to impose a desired behavior.
So let us choose a possible desired open loop to impose through a positive inner feedback (2x2) decoupler of the form following the control scheme:
inverted decoupling.png
The choices , , , , guarantee that , with the desired behavior of the open loop, decoupled and stabilizated.
The resulting control block scheme is:
final decoupler.png
We choose the desired open loop such that it is diagonal, type 1, with all poles i.e. close loop stable and corresponding to two equivalent and indipendent SISO systems.
Its transfer function is:
Ldes = [100/(s*(s+20)), 0 ;0 ,100/(s*(s+20))]; % inizializing the desired open loop transfer function
 
Dd11 = Ldes(1,1)/ Gs_nmp(1,1);
Dd22 = Ldes(2,2)/ Gs_nmp(2,2);
 
Do12 = - (Gs_nmp(1,2)/(Ldes(1,1)));
Do21 = - (Gs_nmp(2,1)/( Ldes(2,2)));
 
Dd = [Dd11, 0;0, Dd22]; % building Dd(s)
Do = [0, Do12; Do21, 0]; % building Do(s)
 
D_ID = Dd*(inv(eye(2)-(Do*Dd))); % building D(s)
L_ID_app = Gs_nmp * D_ID % computing G(s) D(s)
L_ID_app = From input 1 to output... 3.794e07 s^17 + 3.802e09 s^16 + 1.525e11 s^15 + 3.065e12 s^14 + 3.096e13 s^13 + 1.275e14 s^12 + 2.447e13 s^11 + 1.871e12 s^10 + 6.972e10 s^9 + 1.16e09 s^8 - 1.528e06 s^7 - 3.709e05 s^6 - 5055 s^5 - 4.737 s^4 + 0.4714 s^3 + 0.004267 s^2 + 1.206e-05 s 1: ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ 3.794e05 s^19 + 4.561e07 s^18 + 2.286e09 s^17 + 6.116e10 s^16 + 9.227e11 s^15 + 7.466e12 s^14 + 2.574e13 s^13 + 4.913e12 s^12 + 3.75e11 s^11 + 1.396e10 s^10 + 2.32e08 s^9 - 3.094e05 s^8 - 7.423e04 s^7 - 1011 s^6 - 0.9426 s^5 + 0.09432 s^4 + 0.0008536 s^3 + 2.411e-06 s^2 1.546e-06 s^16 + 0.0001238 s^15 + 0.003719 s^14 + 0.04965 s^13 + 0.2498 s^12 + 0.008789 s^11 - 0.003662 s^10 - 0.0004425 s^9 - 2.05e-05 s^8 - 4.694e-07 s^7 - 4.54e-09 s^6 + 3.638e-11 s^5 + 1.464e-12 s^4 + 1.643e-14 s^3 + 8.674e-17 s^2 + 1.813e-19 s 2: --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 1.141e06 s^19 + 1.371e08 s^18 + 6.871e09 s^17 + 1.838e11 s^16 + 2.771e12 s^15 + 2.24e13 s^14 + 7.701e13 s^13 + 1.342e13 s^12 + 9.03e11 s^11 + 2.798e10 s^10 + 3.234e08 s^9 - 3.443e06 s^8 - 1.406e05 s^7 - 1258 s^6 + 4.036 s^5 + 0.1511 s^4 + 0.001031 s^3 + 2.411e-06 s^2 From input 2 to output... -2.375e-07 s^16 - 1.916e-05 s^15 - 0.0005817 s^14 - 0.007965 s^13 - 0.04272 s^12 - 0.02441 s^11 - 0.003906 s^10 - 0.0003052 s^9 - 1.431e-05 s^8 - 3.278e-07 s^7 - 3.725e-09 s^6 - 2.728e-12 s^5 + 4.974e-13 s^4 + 6.439e-15 s^3 + 3.469e-17 s^2 + 7.962e-20 s 1: ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ 3.794e05 s^19 + 4.561e07 s^18 + 2.286e09 s^17 + 6.116e10 s^16 + 9.227e11 s^15 + 7.466e12 s^14 + 2.574e13 s^13 + 4.913e12 s^12 + 3.75e11 s^11 + 1.396e10 s^10 + 2.32e08 s^9 - 3.094e05 s^8 - 7.423e04 s^7 - 1011 s^6 - 0.9426 s^5 + 0.09432 s^4 + 0.0008536 s^3 + 2.411e-06 s^2 1.141e08 s^17 + 1.143e10 s^16 + 4.585e11 s^15 + 9.211e12 s^14 + 9.293e13 s^13 + 3.817e14 s^12 + 6.69e13 s^11 + 4.508e12 s^10 + 1.398e11 s^9 + 1.618e09 s^8 - 1.718e07 s^7 - 7.027e05 s^6 - 6291 s^5 + 20.14 s^4 + 0.7551 s^3 + 0.005153 s^2 + 1.206e-05 s 2: --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 1.141e06 s^19 + 1.371e08 s^18 + 6.871e09 s^17 + 1.838e11 s^16 + 2.771e12 s^15 + 2.24e13 s^14 + 7.701e13 s^13 + 1.342e13 s^12 + 9.03e11 s^11 + 2.798e10 s^10 + 3.234e08 s^9 - 3.443e06 s^8 - 1.406e05 s^7 - 1258 s^6 + 4.036 s^5 + 0.1511 s^4 + 0.001031 s^3 + 2.411e-06 s^2 Continuous-time transfer function. Model Properties
Since we have the usual numerical problems, we substitute zero at the off-diagonal positions.
Let us now compare the resulting open loop with the desired one from a frequency domain point of view.
L_ID= [L_ID_app(1,1), 0;0, L_ID_app(2,2)]; % resultingo L(s) with 0 off-diagonal
 
figure
 
subplot(211)
 
bodemag(Ldes); grid on
Tit = title('Magnitude of Desired Open Loop');
set(Tit,'fontName','Courier','FontSize',10);
xlim([1E-0,1E+2]);
 
subplot(212)
 
bodemag(L_ID);grid on
Tit = title('Magnitude of Obtained Open Loop');
set(Tit,'fontName','Courier','FontSize',10);
xlim([1E-0,1E+2]);
figure
 
subplot(211)
 
phdes = bodeplot(Ldes); grid on
setoptions(phdes,'MagVisible', 'off');
Tit = title('Phase of Desired Open Loop');
set(Tit,'fontName','Courier','FontSize',10);
 
subplot(212)
 
phobt = bodeplot(L_ID); grid on
setoptions(phobt,'MagVisible', 'off');
Tit = title('Phase of Obtained Open Loop');
set(Tit,'fontName','Courier','FontSize',10);
First thing we note is that while the two SISO model on the diagonal of the obtained open loop exhibits the same phase, they differ a bit from the point of view of the magnitude, signal that we are made little numerical errors in computations.
The form of the phases and of the magnitudes are more or less the same as the ones of the desired open loops, but they result a bit anticipated in frequency, exept for the magnitude of the second element on the diagonal, that it exactly the same as the desired one.
Let us evaluate the performances in terms of stability and reference reproduction's capabilities.
T_ID = L_ID(1,1)/(1+L_ID(1,1)); % complementary sensitivity
r_unit3 = [ones(length(t),1)]; % dummy step reference to check if the controller works
 
figure()
 
ycl3 = lsim(T_ID, r_unit3, t);
plot(t, ycl3, 'LineWidth', 2), grid on;
title("Steady state for each decoupled I/O channel");
xlabel("t");
ylabel("y(t)");
legend('y1(t)= y2(t)','FontName','courier','FontSize',12,'Location','SouthEast')
xlim([0 15]);
ylim([-0.5 1.3])
So as expected the inverted decoupler is able to make the closed loop system decoupled and exhibiting a behavior close to the desired, although the overall open loop is not excactly the same as the desired one.
We expect an high control effort peak due to the fact that we are decoupling and stabilizing at the same time.
Cso = D_ID/(eye(2)+ L_ID_app); % control sensitivity function
 
m0 = lsim(Cso, r_unit2, t);
plot(t,m0,'LineWidth',2); grid
title('Control input for r_{i}(t) = step with Inverted Decoupler','FontName','courier','FontSize',12)
xlabel('time','FontName','courier','FontSize',14);
xlim([0 15]);
ylabel('u_{i}(t)','FontName','courier','FontSize',14);
legend('u1(t)',' u2(t)','FontName','courier','FontSize',12,'Location','NorthEast')